Advances in neural engineering now enable modeling of selected pathological aspects of neurodegenerative disease in human neural networks in vitro (Levi et al., 2021; Valderhaug et al., 2021; Estévez-Priego et al., 2023). Of particular interest is how neurodegenerative pathology may manifest and affect network behavior at the microscale and mesoscale, i.e., before it progresses to affect the entire system. The relevant questions are highly challenging or per definition not feasible to address using in vivo models, thus advanced neural engineering models represent a powerful alternative or complementary approach. A fundamental property of such models is that neurons in vitro maintain intrinsic behavior of neurons in the brain by self-organizing over time into networks of increasing structural and functional complexity (Innocenti and Price, 2005; Chiappalone et al., 2006; Giugliano et al., 2006). By combining custom-designed microfluidics devices (MFDs) with microelectrode arrays (MEAs), it is possible to structure multi-nodal human neural networks with controllable connectivity and study both intrinsic network behavior as well as responses to inherent or selectively induced perturbations (van de Wijdeven et al., 2018, 2019; Valderhaug et al., 2019, 2021; Fiskum et al., 2021, 2024; Levi et al., 2021; Heiney et al., 2022; Winter-Hjelm et al., 2023; Weir et al., 2023a,b).

Mutations in the leucine-rich repeat kinase 2 (LRRK2) gene are linked to both late-onset familial and sporadic forms of Parkinson’s disease (PD) (Healy et al., 2008). Although rare, LRRK2 gene mutations have been termed a potential “Rosetta stone” of parkinsonian disorders as all of the major pathologies related to parkinsonism have been observed, in addition to there being end-stage variability, within families carrying the same pathogenic variant (Trinh et al., 2006). The LRRK2 gene is expressed both in the brain and in other tissues throughout the body and is translated into the LRRK2 protein, which has enzymatic kinase activity involved in a range of cellular processes (Trinh et al., 2006). In a recent proteome wide association study of LRRK2 variants from the Parkinsons progression markers initiative cohort (PPMI) (Marek et al., 2018) the LRRK2 gene was directly linked to the regulation of several PD-associated proteins, many of which were subsequently identified as enriched in specific endolysosomal pathways, in microglial cells, and in immune response (Phillips et al., 2023). Among the LRRK2 mutations, the particular G2019S variant represents the most commonly identified cause of late-onset PD, and has been shown to contribute uniquely to both familial and sporadic forms of the disease (Di Fonzo et al., 2005; Gilks et al., 2005; Kachergus et al., 2005; Lesage et al., 2005; Nichols et al., 2005; Ozelius et al., 2006; Trinh et al., 2006; Healy et al., 2008; Moore, 2008; Bouhouche et al., 2017). Moreover, this mutation has been shown to result in differential regulation of a variety of cellular pathways highly relatable to several aspects of PD pathology, among which are axonal guidance, mitochondrial function and calcium homeostasis, cytoskeletal transport, lysosomal function, cell growth, differentiation, and synaptic function (Habig et al., 2008, 2013; Ludtmann et al., 2019; Aleknonytė-Resch et al., 2023; Ohtonen et al., 2023).

We hypothesized that neural networks carrying the Parkinson’s related LRRK2 G2019S mutation would display inherently different micro-and mesoscale behaviors compared to isogenic healthy controls, both during development and in response to induced perturbation. To address this hypothesis, we engineered multi-nodal neural networks from iPSC-derived human neural stem cells with and without the LRRK2 G2019S mutation using a custom-designed MFD (van de Wijdeven et al., 2018, 2019), allowing for investigation of structure–function relationships. Specific microtopographies embedded in the MFD, including axon tunnels and synaptic compartments, allowed for the selective manipulation and targeted investigation of neurites and their organelles, as well as their synaptic elements, i.e., network sites particularly related to early manifestation of neurodegenerative disease pathology. Furthermore, the electrophysiological profile of the multi-nodal neural networks from each group was evaluated in parallel assays by combining the MFDs with an MEA interface.

Monitoring the early network structure and function revealed distinct differences both at the micro-and mesoscale between the LRRK2 G2019S and healthy control neural networks. Specifically, we found that the LRRK2 G2019S neurons self-organized into networks with heightened metabolic demand – exhibiting aberrant morphology and containing neuritic mitochondria moving at higher speeds, compared to the healthy isogenic controls. Inducing a transient, confined overexcitation in targeted nodes of the neural networks further revealed differences in mitochondrial motility and neuritic remodeling between the groups. Together with the indicated difference in total network correlation profiles 24 h post perturbation, these dissimilarities point toward a resilience to change in the LRRK2 G2019S neural networks.

Studies using both in vitro and in vivo models of PD suggest that axonal dysfunction and synaptic alterations represent the earliest detectable signs of the disease (MacLeod et al., 2006; Cheng et al., 2010; Dagda et al., 2014) and that initiation of pathology at the axon terminals might signify the start of the retrograde degeneration of the neurons (Tagliaferro and Burke, 2016; Sheehan and Yue, 2018). Although extensively reported on Smith et al. (2005), West et al. (2005), MacLeod et al. (2006), Smith et al. (2006), Plowey et al. (2008), Gillardon (2009), Dachsel et al. (2010), Meixner et al. (2011), Cherra et al. (2013), Habig et al. (2013), and Sepulveda et al. (2013), the specific involvement of the LRRK2 gene in neurite process morphology remains obscure. We found the LRRK2 G2019S neural networks to have a much greater volume of neurites within the axonal tunnels and synaptic compartments, as well as random outgrowth profiles compared to the healthy isogenic control. This result was corroborated through several different imaging approaches, i.e., immunocytochemistry, bright field microscopy, Calcein-AM labeled cells, as well as live imaging of TMRM labeled mitochondria. In contrast to our findings, the specific G2019S PD-associated LRRK2 mutation has generally been found to result in a progressive reduction of neurite length and branching (West et al., 2005; Greggio et al., 2006; MacLeod et al., 2006; Smith et al., 2006; Plowey et al., 2008; Chan et al., 2011; Nguyen et al., 2011; Winner et al., 2011; Sanchez-Danes et al., 2012; Cherra et al., 2013; Reinhardt et al., 2013; Dagda et al., 2014; Qing et al., 2017), due to an increase in kinase activity (with one exception demonstrating non-impaired neurite morphology (Dachsel et al., 2010)). Moreover, knockdown and knock-out models resulting in LRRK2 deficiencies present with a progressive increase in neurite length in some studies (MacLeod et al., 2006; Dachsel et al., 2010; Winner et al., 2011; Habig et al., 2013; Sepulveda et al., 2013), while others find the opposite (Gillardon, 2009; Meixner et al., 2011).

The morphology of neurites and plasticity of synapses is affected by mitochondrial function as the availability of mitochondria is both essential and limiting for the support and maintenance of these structures (Li et al., 2004). A loss or impairment in dendritic mitochondria is a central feature of PD pathogenesis (Bose and Beal, 2016; Verma et al., 2017; Singh et al., 2019), and other studies on the particular LRRK2 G2019S variant have indicated compromised mitochondrial function (Cooper et al., 2012; Sanders et al., 2014; Hsieh et al., 2016; Schwab et al., 2017; Bose and Beal, 2019; Ludtmann et al., 2019). Here, the mitochondrial dynamics analysis showed that the LRRK2 G2019S neural networks contained about twice the ratio of motile mitochondria (mean 0.139 vs. 0.06), moving at about twice the speed (mean 1.296 μm/s vs. 0.639 μm/s), and with clear directional tendencies, compared to the mitochondria of the healthy isogenic control neural networks. Our baseline measurements thus point toward an overall higher metabolic cost of LRRK2 G2019S neural networks compared to healthy isogenic controls. These measurements are in line with the results from two other comparable studies, which also found an increase in mitochondrial motility in iPSC-derived neurons carrying the G2019S LRRK2 mutation compared to controls (Cooper et al., 2012; Schwab et al., 2017). Although their study was performed after 4–5 weeks of differentiation, Cooper et al. also reported an increase in mitochondrial velocity, with the LRRK2 mitochondria moving at about 0.7 μm/s versus 0.4 μm/s for the controls (Cooper et al., 2012). Moreover, in a study by Othonen and colleagues, RNA sequencing of microglia-like-cells (iMLCs) derived from iPSCs carrying the LRRK2-G2019S mutation showed alterations in pathways related to oxidative stress, which resulted in deficits in mitochondrial respiration (Ohtonen et al., 2023). In further support of our findings, these authors also report a higher oxygen consumption rate in LRRK2-G2019S carrying iMLCs compared to isogenic controls, resulting in higher respiration at the basal state.

Interestingly, previous studies investigating neuronal axonal transport kinetics, both in vitro and in vivo (Lewis et al., 2016), as well as in networks structured using microfluidic devices (Moutaux et al., 2018), have found mitochondrial dynamics to decrease with neural network maturation and development, and mitochondrial immobilization and stability to be important hallmarks of mature axons. Thus, the in vitro neural network traits displayed by the LRRK2 G2019S networks in our study are more prominent at very early time points in general network development. This might further suggest an impairment in network development, a notion supported by our observation of differential growth cone profile at very early stages of LRRK2 G2019S neuron differentiation compared to controls. Altered neural growth cone morphology and number has been found by others after knocking down LRRK2 (Habig et al., 2013). Importantly, the growth cone has pronounced influence on mitochondrial transportation, as nerve growth factor (NGF) signaling from immature axons regulates mitochondrial motility (Morris and Hollenbeck, 1993; Chada and Hollenbeck, 2003; Chada and Hollenbeck, 2004; Saxton and Hollenbeck, 2012). Furthermore, impaired axonal growth and guidance would explain what is observed as aberrant (and consequently inefficient) network wiring in the LRRK2 G2019S networks, with large numbers of neurites crossing perpendicular to the axonal tunnels in the synaptic compartment. Based on the general structure–function relationships of neural networks (Valderhaug, 2020), resulting aberrant wiring might again affect neurotransmission efficacy, causing less efficient signal propagation within the network, again affecting the overall energy status of the neuron and consequently the motility and distribution of mitochondria (Saxton and Hollenbeck, 2012).

Taken together, the baseline results indicate that the LRRK2 G2019S neurons self-organize into multi-nodal neural networks with inefficient wiring (aberrant neurite morphology) as well as increased mitochondrial dynamics, relative to the healthy isogenic control neural networks. Thus, both micro-and mesoscale baseline features point toward an overall heightened metabolic cost associated with the LRRK2 G2019S mutation, suggesting a network phenotype more vulnerable to perturbation.

The multiple-hit hypothesis of PD suggests that the development of the disease results from a combination of factors, rather than a single cause. This hypothesis aligns with the understanding that PD is a complex neurodegenerative disorder with a multifactorial etiology, encompassing genetic predispositions, environmental exposure, and their interactions may contribute to the onset and development of the disease (Carvey et al., 2006; Sulzer, 2007; Patrick et al., 2019). Based on our baseline results and the variability in disease progression in patients with LRRK2 G2019S associated PD, it is reasonable to assume that some phenotypic expression of the mutation may become apparent only following a significant or stressful challenge (Benson et al., 2018). Following the transient perturbation (KA) in our study, the healthy isogenic control neural networks revealed greater responses in almost all measures: a significant immediate reduction in number of neuritic mitochondria, as well as a more prominent modification in neurite morphology and synaptic remodeling, and the only measured decrease in total network correlation, 24 h post stimulation, relative to the LRRK2 neural networks. The mitochondrial motility was the only measure in which the LRRK2 neural networks displayed a greater response, with a significant reduction in mitochondrial motility immediately following perturbation, while the control neural networks showed no change.

The highly significant difference found in neurite morphology in the control neural networks with a retraction/reduction of boutons observed in the synaptic chambers 24 h post KA stimulation, suggests neurite remodeling in response to the transient perturbation. At the same time, the number of synapses (co-occurrence of pre- and post-synaptic markers) in the control neural networks was not significantly altered, but the size of their overlapping area was, with much larger synaptic areas measured 24 h after overexcitation. This rapid activity-dependent alteration in dendritic spine morphology is likely an expression of a regular mechanism for converting short-term synaptic activity to long term lasting changes in connectivity and function. In line with our baseline measurements, the size range of the postsynaptic density is usually within 0.2–0.5 μm, and can be localized to both spiny and non-spiny structures. Furthermore, this area contains both the kainate and AMPA receptors (Sheng, 2001), i.e., glutamatergic receptors targeted by our stimulation. During synaptic plasticity, the PSD increases in size in response to potentiation events, and the glutamate receptors contained within can be modulated by neural activity on a timescale from minutes to weeks (Sheng, 2001). In contrast to the observed structural changes in the control neural networks, no significant alteration in synaptic number or size was found in the synaptic compartment of the LRRK2 mutated neural networks 24 h after KA stimulation, suggesting impaired synaptic plasticity.

Importantly, this goes hand in hand with our electrophysiology results, where the KA-stimulated isogenic control neural networks were the only ones to display a relative decrease in total network correlation at the 24 h post timepoint as well, while all other conditions displayed a relative increase. Measures of MFR were comparable at all analyzed timepoints. This is in line with the findings of another highly relevant study, where the average neuronal activity per se was found not to be sufficiently informative to reveal disease-related alterations (Carola et al., 2021). As the transient excitatory event indiscriminately excites related connections that are both functional and non-functional, and likely produces activity-dependent synaptic modifications at both local and distal network sites (Bi and Poo, 2001), a decrease in total correlation could be expected after a “healthy” plastic network response to such an event.

Furthermore, other studies have found increased vulnerability to oxidative stress, higher levels of mtDNA damage, and impaired mitochondrial movement as a result of the G2019S mutation, indicating compromised mitochondrial function (Cooper et al., 2012; Sanders et al., 2014; Hsieh et al., 2016; Schwab et al., 2017; Bose and Beal, 2019). Our baseline measurements indicated mitochondrial activity of greater metabolic cost for the LRRK2 neural networks, which in turn would make these networks more vulnerable to perturbation. As already noted, alterations in mitochondrial function in turn affect the plasticity of synapses and morphology of neurites (Li et al., 2004). The energy status of the neuron greatly affects the motility and distribution of mitochondria (Saxton and Hollenbeck, 2012), and a relocation of the mitochondria toward either the soma or the presynaptic terminal during or following an overexcitation event is to be expected as this corresponds to the locations of greatest metabolic demand at the time. This fits with our investigations, where a reduction in number of mitochondria was found in the neurites of the control neural networks immediately following the induced excitatory event. The motility ratio of these mitochondria in turn remained unchanged, likely reflecting a regular response with mitochondria being recruited to other areas with greater metabolic demand. The LRRK2 neural networks on the other hand showed the opposite result, with no change in number of mitochondria and a significant decrease in mitochondrial motility immediately after the excitatory event. This fits well with our finding of substantial neurite and synaptic remodeling in the control neural networks only, as the mitochondria of the LRRK2 neural networks seemingly become less mobile in response to the perturbation rather than being efficiently recruited for energy metabolism and remodeling at the stimulated synapses. This also aligns well two recent in vivo studies, where the LRRK2 G2019S mutation was linked to impaired experience-dependent plasticity (Matikainen-Ankney et al., 2018), abnormal synaptic changes and a lack of adaptive change in intrinsic excitability (Guevara et al., 2020), in response to social stress.

The reduced motility and unaltered number of active mitochondria observed in the neurites of the LRRK2 neural networks could be a consequence of the aberrant structural network wiring, with altered signal propagation leading to impaired recruitment of mitochondria, and/or it could be due to impaired mitochondrial transport and/or function. Nevertheless, together, the baseline elevated metabolic cost and the lack of both structural and functional network alterations following perturbation, produces a network seemingly resilient to change, associated with the LRRK2 G2019S mutation. In the long term, this inefficient and energy-demanding phenotype renders the network more vulnerable.

Although graph-theory based analyses were beyond the scope of the current study, our findings are highly relevant in the context of other studies where this type of analysis was applied to either electrophysiology or fluorescent calcium imaging data. Exaggerated small-worldness, in particular, is a network structure associated with high metabolic cost, and inefficient, noisy information transfer between network regions (Bassett and Bullmore, 2017). Prominent small-world organization has been identified to underly the metabolic patterns and pathological alterations in brain network structure and function of PD-patients using neuroimaging data (Niethammer and Eidelberg, 2012; Ko et al., 2018; Schindlbeck and Eidelberg, 2018; Schindlbeck et al., 2019). Additionally, recent findings by our group using in vitro electrophysiological recordings from human neural networks with inherent neurodegenerative pathology show that small-world propensity associates with increased metabolic demands (Fiskum et al., 2024). This renders neural networks with predisposing genetic mutations particularly vulnerable to neurodegenerative processes, which can be precipitated by transient external perturbation. Furthermore, when compared to GBA variants of PD, the phenotypically slower progression of disease at early stages in LRRK2-PD can be explained by the gain of functional connections found in central network cores (linking the cerebellum and putamen), which could compensate and provide robustness against early symptomatic expression in prodromal stages (Schindlbeck et al., 2019). In line with this, in a relevant in vitro study by Carola and colleagues (Carola et al., 2021), information theory analysis was applied to fluorescence calcium imaging data. It showed that LRRK2 PD neural networks to have large, strongly linked functional communities, with lower average connectivity, ultimately producing a tendency toward greater synchrony. Thus, the overall network resilience to change in response to perturbation displayed by the LRRK2-G2019S neural networks in our study, compared to the healthy controls, reflects both underlying impairments and compensatory mechanisms.

In summary, our study points toward an overall heightened metabolic demand, with aberrant morphology and mitochondrial dynamics, related to the G2019S LRRK2 mutation. Furthermore, these alterations seem to cause a network resilience to perturbation at early time points, as transient excitation revealed a lack of neurite remodeling- and synaptic plasticity response in LRRK2 G2019S mutated neural network compared to healthy isogenic controls. Whether such a resilience is adaptive or maladaptive, compensatory or pathological might depend on the experimental timeframe, and might cause different functional outcomes depending on PD disease stage. Thus, advanced multidisciplinary approaches such as ours, where several relevant methodologies are combined, are needed if we are to progress the field. The ability to recapitulate the relevant structural and functional dynamics of these networks at the microscale and mesoscale level opens up for entirely new avenues in modeling PD, and by the same token, other neurodegenerative diseases, including amyotrophic lateral sclerosis and Alzheimer’s disease. As such, the advanced modeling approach and new findings presented in this study are highly relevant in the quest for elucidating underlying disease mechanisms. Equally importantly, the relevant insights are fundamental for formulating new hypotheses, as well as identifying time window(s) and mode(s) of therapeutic intervention, with a view to clinical translation.

Interestingly, in a recent seminal study using alpha-synuclein seed amplification assay (SAAs) to diagnose PD with a staggering 93% accuracy, the LRRK2 PD subgroup had one of the lowest proportions of positive SAA test (67.5%) (Siderowf et al., 2023). This finding echoes the 1/3 frequency of LRRK2 PD individuals reported to lack alpha-synucleinopathy (Lewy pathology) in post-mortem studies (Kalia et al., 2015; Schneider and Alcalay, 2017), and represents a peculiarity in disease etiology which warrants further elucidation. In this context, we have shown in a previous study that this approach can be used to investigate the early functional changes associated with alpha-synuclein proteinopathy in engineered human neural networks (Valderhaug et al., 2021). It is thus highly feasible and relevant to combine such advanced multidisciplinary approaches to investigate the micro-and mesoscale aspects of LRRK2 human neural networks with a focus on the effect and development of alpha-synucleinopathy in the future. Furthermore, analogs of the validated PD network biomarker – PDRP (Parkinsons disease related pattern) (Niethammer and Eidelberg, 2012) and PDCP (Parkinsons disease cognitive pattern) (Schindlbeck and Eidelberg, 2018)– could be applied to the electrophysiological data using graph theory, and utilized to investigate disease specific network patterns, such as small-worldness (Bassett and Bullmore, 2017; Ko et al., 2018; Schindlbeck et al., 2019; Valderhaug et al., 2019; Valderhaug, 2020; Fiskum et al., 2024), in advanced cellular models based on engineered neural networks. Furthermore, fluorescence imaging can be applied to enable to visualize neuronal activity for functional characterization in such disease models. Such an approach was successfully applied in the aforementioned study by Carola and colleagues (Carola et al., 2021). Using information theory to analyze calcium activity data, revealed distinct activity patterns separating iPSC derived neuronal cultures from those of LRRK2 PD patients from isogenic controls and healthy donors. Furthermore, combination of calcium imaging with MEA-based electrophysiology can in future studies provide a more comprehensive profile of network dynamics at high spatial and temporal resolution, respectively. Such modeling approaches may thus provide a means to pinpoint structural and functional network sources that might be causally related to phenotypic variation, for instance in treatment response.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

The human cell lines utilised in this study are commercially available and their use is regulated by a Material Transfer Agreement between the supplier and NTNU. All cell work was carried out in GMO/GMM approved Biofacility level 2 laboratories. Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements.

VV: Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. OR: Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. RW: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing. KH: Formal analysis, Methodology, Writing – original draft, Writing – review & editing. SN: Supervision, Writing – original draft, Writing – review & editing. AS: Conceptualization, Project administration, Supervision, Writing – original draft, Writing – review & editing. IS: Conceptualization, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing.