Developmental effect of RASopathy mutations on neuronal network activity on a chip

RASopathies are a group of genetic disorders caused by mutations in genes encoding components and regulators of the RAS/MAPK signaling pathway, resulting in overactivation of signaling. RASopathy patients exhibit distinctive facial features, cardiopathies, growth and skeletal abnormalities, and varying degrees of neurocognitive impairments including neurodevelopmental delay, intellectual disabilities, or attention deficits. At present, it is unclear how RASopathy mutations cause neurocognitive impairment and what their neuron-specific cellular and network phenotypes are. Here, we investigated the effect of RASopathy mutations on the establishment and functional maturation of neuronal networks. We isolated cortical neurons from RASopathy mouse models, cultured them on multielectrode arrays and performed longitudinal recordings of spontaneous activity in developing networks as well as recordings of evoked responses in mature neurons. To facilitate the analysis of large and complex data sets resulting from long-term multielectrode recordings, we developed MATLAB-based tools for data processing, analysis, and statistical evaluation. Longitudinal analysis of spontaneous network activity revealed a convergent developmental phenotype in neurons carrying the gain-of-function Noonan syndrome-related mutations Ptpn11^D61Y and Kras^V14l. The phenotype was more pronounced at the earlier time points and faded out over time, suggesting the emergence of compensatory mechanisms during network maturation. Nevertheless, persistent differences in excitatory/inhibitory balance and network excitability were observed in mature networks. This study improves the understanding of the complex relationship between genetic mutations and clinical manifestations in RASopathies by adding insights into functional network processes as an additional piece of the puzzle.

Introduction

RASopathies cover a group of genetic disorders associated with a characteristic pattern of anomalies, including distinctive facial features, cardiopathies, growth and skeletal abnormalities and various degrees of developmental delay (Zenker and Kutsche, 2016). They are typically caused by germline mutations in genes encoding components or regulators of the RAS/mitogen-activated protein kinase (MAPK) signaling pathway. Common is an overactivation in signaling, either by activating mutations in the pathway components or by loss of function in negative modulators (Tartaglia et al., 2010; Zenker and Kutsche, 2016). The RAS/MAPK signaling pathway has a central role in the cellular response to extracellular stimuli (Molina and Adjei, 2006; Yoon and Seger, 2006). In neurons, it is indispensable for proper differentiation, synapse formation and plasticity and is therefore critical for brain development, learning, memory and cognition (Samuels et al., 2009). Indeed, most RASopathies are associated with neurocognitive impairments and variable degrees of intellectual disability (Cesarini et al., 2009; Tajan et al., 2018; Kang and Lee, 2019). While patients with cardio-facio-cutaneous syndrome (CFC), Costello syndrome (CS) and SYNGAP1-related encephalopathy show moderate to severe intellectual disabilities and neurodevelopmental delay, epilepsy and autism spectrum disorder (Roberts et al., 2006; Armour and Allanson, 2008; Schulz et al., 2008; Adviento et al., 2014; Pierpont et al., 2022), patients with Noonan syndrome (NS) and neurofibromatosis type 1 are generally less severely affected. Common diagnoses are mild neurodevelopmental delay, impaired executive functioning, learning difficulties and attention deficit hyperactivity disorder symptomatology (Sharland et al., 1992; Mautner et al., 2002; Krab et al., 2008; Payne et al., 2011; Alfieri et al., 2014; Pierpont et al., 2015; Kim and Baek, 2019).

Interestingly, some of the neurocognitive features diminish with increasing age of the patients, while different features emerge. For example, children with NS have shown more extensive cognitive problems concomitant with learning difficulties and intelligence impairment compared to adults (Roelofs et al., 2016).

In animal experiments, the development of neural networks has not yet been studied in a systematic way. While most studies demonstrated electrophysiological and behavioral phenotypes in mice carrying different RASopathy mutations, most of them did not assess the respective phenotypes along the animal lifespan (Costa et al., 2002; Cui et al., 2008; Lee et al., 2014; Altmuller et al., 2017). An exception is a recent study that used mice in which the Noonan syndrome-related Kras^G12V mutation was expressed exclusively in neurons. In the hippocampi of these animals, upregulation of RAS/MAPK signaling was observed during the early phases of postnatal development but not in the adult state (Papale et al., 2017). Interestingly, the normalization of RAS/MAPK pathway activity was concomitant with enhanced GABAergic synaptogenesis, suggesting that maladaptive development of the neural network may underlie changes in neurocognitive effects of RASopathy observed in patients across the lifespan. Importantly, the switch in the cellular mechanism underlying network dysfunction might also have consequences for the selection of an effective treatment strategy. In line with that, Papale et al. (2017) were able to improve the cognitive phenotypes in Kras^G12V mice by treatment with RAS/MAPK inhibitors applied in juvenile animals. Strikingly, the same treatment was not effective if applied in older animals. Here, in line with increased GABAergic synaptogenesis, inhibition of GABAergic neurotransmission showed positive effects.

In this study, we test the hypothesis that pathological activation of RAS/MAPK signaling may induce maladaptive processes leading to compensatory effects in developing neuronal networks. In vivo investigation by real-time monitoring of neuronal activity and functional connectivity is technically challenging and not suitable for large-scale testing of therapeutic compounds (Bang et al., 2021). As an alternative opportunity to study functional connectivity and network parameters, we investigated neuronal cell cultures derived from RASopathy mouse models plated on multiwell MEA (mwMEA) plates in a longitudinal study. We used previously characterized mouse models of Noonan syndrome that bear gain-of-function mutations in Ptpn11^D61Y (Altmuller et al., 2017) or in Kras^V14l (Schubbert et al., 2006; Gremer et al., 2011; Hernandez-Porras et al., 2014). To perform an extensive functional analysis of spontaneous network activity and electrically evoked activity, we used an in-house developed analysis tool written in MATLAB and routines for the application of principal component analysis (PCA) to reduce dimensionality in the data set for better data visualization and statistical evaluation. We observed significant effects of RASopathy mutation on network activity at early developmental stages that were compensated during network maturation. However, the mature RASopathy networks with apparently normal network activity still exhibited differences in excitatory/inhibitory balance and network excitability.

Materials and methods

Animals

For this study, conditional Ptpn11^{D61Yfloxed/wt} mice (B6.129S6-Ptpn11tm1Toa/Mmjax, Jackson Laboratories, RRID# MMRRC_032103-JAX) (Chan et al., 2009) and Kras^{V14lfloxed/wt} mice (Hernandez-Porras et al., 2014) were used. In both strains, a LOX-STOP-LOX cassette is cloned before the exon containing the RASopathy mutation. Constitutive expression of these mutations leads to severe somatic phenotypes (Araki et al, 2009; Chan et al., 2009; Hernandez-Porras et al., 2014). To reduce the burden of experimental animals, we generated animals with expression of RASopathy mutations restricted to neocortex. To this end heterozygous mice bearing the mutation were crossed with mice homozygous for Emx1-cre alleles (B6.129S2-Emx1tm1(cre)Krj/J animals, Jackson Laboratories, RRID# IMSR_JAX:005628, genotyping according to provider’s recommendations). Litters from this cross were all heterozygous for the Emx1-cre allele, controls were homozygous wt, while mutants were heterozygous for the Ptpn11^D61Yfloxed or Kras^V14lfloxed allele. All mouse strains were bred on a C57BL/6N background for more than 5 generations in our laboratory. Mice were kept at 22 ± 2°C on a 12 h light-dark cycle with food and water ad libitum. Breeding of animals and experiments using animal material were performed in accordance with the local animal welfare officer (FAU: TS12/2016) and in accordance with the European Directive 2010/63/EU. Animal records were kept by Python-based Relational Animal Tracking Software (PYRAT-Scionics Computer Innovation GmbH, Dresden, Germany).

Animals were genotyped using published protocols. Specifically, the EMX1-cre genotyping protocol provided by the Jackson Laboratory was followed. To detect the Ptpn11^D61Y allele, polymerase chain reaction (PCR) using TGGAGCTGTTACCCACATCA and GCACAGTTCAGCGGGTACTT primers followed by a melting point analysis using the High Resolution Melting and Gene scanning application on the LightCycler 480 (Roche Diagnostics) was performed according to Altmuller et al. (2017). To determine the genotype in Kras^V14l, PCR was performed using AGG GTA GGT GTT GGG ATA GC, CTC AGT CAT TTT CAG CAG GC, CTG CTC TTT ACT GAA GGC TC primers using S7 Fusion High-Fidelity DNA polymerase (Biozym, Cat#MD-S7-100) according to Hernandez-Porras et al. (2014). The protocol was generated to discriminate between 403 bp (wild type) and 621 bp (knock-in mutation) alleles. Primers were purchased from Eurofins, Ebersberg, Germany.

Dissociated cell culture and plating on multiwell MEA plates or 18 mm coverslips

Chemical reagents are listed in Supplementary Table S1. The day before dissection, wells of 48-well mwMEA plates (Cytoview Cat#MEA 48, M768-tMEA_48B, Axion Biosystems) and 18 mm Menzel glass coverslips (#6311342, VWR International, Radnor, United States) were coated with poly-L-lysine (0.5 mg/mL) followed by overnight incubation at 37°C. The next day, the wells and coverslips were washed three times with sterile double-deionized water and left in HBSS−/−, which was removed directly before cell seeding. Prior to seeding of cells on 18 mm coverslips, they were coated with 100 μL of Neurobasal^™-A media (NBA) supplemented with 0.2 mM Glutamax, 2% (v:v) B27, 0.1 M Antibiotic-Antimycotic (named NBA mix) and 10% fetal calf serum (FCS), as previously described (Anni et al., 2021), for better attachment of neurons.

Dissociated neuronal cultures were prepared from individual brains of newborn mice (P0–P1) as previously described in Anni et al. (2021). Neonates of both sexes were used for culture preparation. Briefly, Ptpn11^{D61Yfloxed/wt} (here termed Ptpn11^D61Y) or Kras^{V14lfloxed/wt} (here termed Kras^V14l) mice and their siblings (here termed control) were decapitated, the brains were removed, and the forebrains were separated and freed of meninges. Dissociation of cells took place chemically in Papain-protease-DNase mix (HBSS−/−, 0.01% Papain reconstituted in EBSS, 0.01% (w:v) DNase I, 0.1% (w:v) Dispase II) at 37°C for 10 min followed by mechanical dissociation. This process was repeated 2 times. The cell suspension was passed through a 70 μm cell strainer. Cells were spun down and washed with NBA mix three times and then resuspended in NBA. Cells were counted and diluted to plate 80,000 cells in a volume of 50 μL onto the coated wells of the mwMEA plate. The cells were incubated at 37°C in a 5% CO₂ atmosphere for 1 h and allowed to settle, followed by the addition of 250 μL of NBA mix for further culturing. Similarly, 200,000 cells in a volume of 100 μL were carefully seeded onto coverslips for 1 h immediately after removing the NBA mix containing FCS droplets. They were then transferred to 12-well plates containing 1 mL of NBA mix per well and maintained at 37°C in 5% CO₂ until maturation. Genotyping was performed post hoc from tail biopsies obtained from euthanized pups.

Immunocytochemistry and image acquisition and analysis

Neuronal cultures were prepared independently form 4 animals per genotype and grown on coverslips for 21 days. From each animal 2 independent coverslips were stained for each antibody combination. Coverslips with neurons from both genotypes were processed in parallel with identical antibodies, solutions, and other chemicals as follows: cells were fixed with 4% (w:v) paraformaldehyde in PBS for 4 min at room temperature (RT), washed in PBS and blocked and permeabilized in tandem with PBS solution containing 10% (v:v) FCS, 0.1% (v:v) glycine and 0.3% (v:v) Triton X-100 for 45 min at RT (Supplementary Table S1). Coverslips were incubated with primary antibodies at 4 degrees overnight and after washing with PBS with secondary antibodies at RT for 1 h. All antibodies were applied in PBS containing 3% (v:v) FCS. Coverslips were mounted on glass slides with Fluoroshield. All antibodies used are listed in the Supplementary Table S1. Immunofluorescence staining of vesicular GABA transporter (VGAT) and vesicular glutamate transporter (VGLUT1) was performed as double staining in combination with integral synaptic vesicle protein synaptophysin. Image acquisition was performed by an epifluorescence microscope (Nikon Eclipse Ti, Nikon Corporation) equipped with an iXon EM+ 885 EMCCD Andor camera (Andor Technology) using a 60X/NA 1.2 objective (Plan APO VC Nikon CFI, Nikon Corporation) and controlled by VisiView software (Visitron Systen GmbH). Identical camera and illumination settings were applied during image acquisition for coverslips of all experimental conditions imaged on the same day. Three to four images were obtained from each cover slip. Data were stored as 16-bit images.

Sixteen-bit images were analyzed using the in-house MATLAB routine SynEval, that includes an iterative, water-shed-based segmentation algorithm and includes preprocessing by background subtraction and normalization (Guhathakurta et al., 2022). SynEval applies filtering, dilation and flood filling to get a rough segmentation within which the threshold is defined as the median of the lowest 1% fluorescence intensity. Here, the pictures were segmented according to immunofluorescence staining for synaptophysin to obtain synaptic regions of interest (ROIs). Subsequently, ROI coordinates were transferred to the channel in which VGAT or VGLUT immunofluorescence was recorded, and then the ROI were identified by being positive for VGAT or VGLUT and the mean immunofluorescence intensity (IFI) within the transferred ROI was obtained. Mean IFI represented the average intensity value derived from pixels within the ROI. Post-processing calculations yielded the IFI of co-labeled puncta and their fraction of the total number of synaptic puncta detected.

MEA recordings

The wells in 48-well MEA plates were equipped with 16 PEDOT electrodes arranged in a 4 × 4 grid with an electrode spacing of 340 μm, electrode diameter of 50 μm and recording area of 1.1 mm x 1.1 mm. Recordings were performed using a Maestro multiwell MEA recorder (Axion Biosystems) at 37°C in a 5% CO₂:95% air atmosphere. Before recording, neuronal cultures were visually checked for vitality under a light microscope (Supplementary Figure S1A). Prior to each recording during the time series experiments, the MEA plate was kept in the MEA recorder for at least 5 min to allow the neuronal cultures to adjust to ambient conditions and to avoid bias by mechanical disturbances. To acquire data and for spike detection, we used Axion Integrated Studio (AxIS) software (2.4.2.13). The heatmap visualization in the software enabled a first control of neuronal activity online (Supplementary Figure S1B). Voltage potentials were recorded across all electrodes with a sampling rate of 12.5 kHz (Supplementary Figure S1C). Recordings were subdivided into 5 min time units to avoid a large amount of data at once (raw file). The raw signal was filtered by a Butterworth filter with a 200 Hz to 2.5 Hz bandpass filter. The filtered signal was detected for neuronal events when the signal exceeded a threshold denoted as spikes (cf. Ch. data analysis, Supplementary Figure S1D). Specifically, an adaptive threshold was applied within a moving window around each spike with 0.84 ms pre-spike and 2.16 ms post-spike range. In this window, the threshold was adaptively set as a multiple of the median noise level [5.5 of standard deviation (SD)] on an electrode (Supplementary Figure S1C). These were stored in a spike list (.csv file) containing information about the time and amplitude of detected spikes chronologically.

Time series of spontaneous network activity

Multiple recordings of spontaneous activity in neuronal networks plated on the same mwMEA plate were longitudinally and noninvasively conducted in the time range 10 days in vitro (DIV) to 33 DIV. On a measuring day, spontaneous activity was recorded for 20 min in 5 min time units.

Disinhibition-induced stimulation experiments

The application of bicuculline (bic) to Ptpn11^D61Y neuronal cultures was performed on 33 DIV. The baseline was recorded for 20 min in 5 min time units, and then the plate was removed to add bic (volume: 10 μL, final concentration 10 μM, diluted in H₂O from a 10 mM stock solution in DMSO). After treatment, the plate was placed back in the MEA recorder and kept there for an additional 10 min without recording, followed by 50 min post bic recordings in 5 min time units (Supplementary Figure S2). Values were normalized to the baseline (100%) and presented as percentages.

Electrical stimulation experiments

Electrical stimulations were applied on mature networks (three independent experiments on DIV29, DIV30 and DIV39). To provide electrical test stimulation (STIM) and tetanus stimulation (TET), individual electrodes from the MEA were triggered to deliver electrical pulses (Supplementary Figure S3A). Protocols for STIM and TET were modified from Chiappalone et al. (2008) (Supplementary Figure S3C). The STIM protocol consists of a train of pulses (0.2 Hz). The pulses were biphasic from +750 mV to −750 mV lasting for 250 μs per phase. During STIM, six electrodes per array (stimulation electrodes e1–6, as shown in Supplementary Figure S3A) were triggered sequentially for 5 min, also denoted as sessions 1–6 (Supplementary Figure S3D). Prior to each STIM, which consequently lasted 30 min in total, spontaneous activity was recorded initially for 20 min (prior STIM1) and subsequently for 10 min (prior STIM2 and STIM3) (Supplementary Figure S3D). TET was delivered from one electrode and consisted of 20 bursts at 0.2 Hz, each containing 11 pulses at 20 Hz (Supplementary Figure S3C).

Analysis of MEA data

Processing MEA-based data

For further data analysis, customized software packages coded in MATLAB were used to process the spike lists (.csv) that were output by AxIS software. Each analysis was run in batch mode to save time and avoid human bias. The packages enable an analysis for each experimental setting, including spontaneous network activity in time series recordings plus network disinhibition experiments as well as electrically evoked activity. We extended the spontaneous network activity analysis for time series experiments with MATLAB programmed functions performing PCA commonly applied in machine learning to reduce dimensionality in data sets.

Spontaneous network activity

The software package for analyzing spontaneous network activity in longitudinal experiments is equipped with a graphical user interface (GUI) (Supplementary Figure S4) to enable the user-friendly input of experiment-related relevant data (cf. scheme in Supplementary Figure S5A). In the first step, the user selects the spike list files (.csv) that had been output by Axion software after each measurement. The minimum spike rate to determine active electrodes (herein 0.1 Hz) and the minimum number of active electrodes (herein, unless otherwise stated: 10) can be specified in the GUI.

Spike calculation

Spikes, which define neuronal firing events, are detected when the filtered voltage signal crosses a threshold in the continuous data stream. The mean firing rate (MFR) is defined as the total number of spikes divided by recording time, MFR = n_spikes/s. The weighted MFR is calculated as wMFR = (e_max/e_act)·MFR. Herein, the maximum number of active electrodes (e_max) was determined at the culture age of maturity (DIV 21–24). e_act refers to the currently active number of active electrodes. An electrode was defined as active when its spiking rate was at least 0.1 Hz, and then it was also described as a contributing channel. For disinhibition experiments, e_max was determined as the number of active electrodes during the baseline measurement prior to treatment. Analogously, wtMFR [mean firing rate weighted to total (t) number of electrodes] indicates the same calculation, but e_max refers to the total number of electrodes (herein 16 electrodes). The interpike interval is the mean time between spikes in s.

Burst detection and calculation

A burst was defined as neuronal activity that appears in high-frequency spiking followed by a period of quiescence on a single electrode level. Herein, burst detection was performed along an interval interspike threshold detecting bursts as spiking events of no longer than 100 ms intervals between 5 consecutive spikes. The burst duration is calculated as the time from the first spike to the last spike within a burst. The mean bursting rate (MBR) is the total number of bursts within one well divided by recording time. The value is weighted analogously to wMFR (wMBR = (e_act/e_max)·MBR).

Network analysis

The analysis of synchronized spiking in network bursts (NB) and correlated spiking measured as spike time tiling coefficient requires the transformation of the spike list to a time stamp matrix, a binary array of sample points indicating spikes as ones and sample points without a spike as zeros. With a sampling rate of 12.5 kHz, one sample point correlates to 0.08 ms. The assessment of network-wide coordinated activity requires the detection of spike clusters across multiple electrodes. Based on the time stamp matrix, a MATLAB-coded routine detected NBs across the array as spike clusters containing a minimum number of 50 spikes with a maximum of 100 ms interspike interval and with at least 5 contributing channels. Then, several NB-describing parameters, such as NB duration, number of NBs per well, and number of spikes per NB, were calculated. The NB duration was calculated as the time from the first spike to the last spike in an NB. To investigate changes in the network dynamics, the correlation of spiking between each pair of electrodes was analyzed by the spike time tiling coefficient (STTC) (Cutts and Eglen, 2014), which reveals the synchronization between pairs of electrodes by returning the unbiased correlation in firing. It is calculated as

$\mathrm{STTC}=\frac{1}{2}\left(\frac{P_{\mathrm{A}}-T_{\mathrm{B}}}{1-P_{\mathrm{A}}{T}_{\mathrm{B}}}+\frac{P_{\mathrm{B}}-T_{\mathrm{A}}}{1-P_{\mathrm{B}}{T}_{\mathrm{A}}}\right)$

where A and B represent two spike trains. P_A (P_B) is the proportion of spikes from A (B) that occurs within the time window ±Δt (100 ms) around any spike from B (A). T_A (T_B) is the proportion of recording time that lies within ±Δt of any spike from A (B). The pairwise computation of STTC between electrodes resulted in a correlation matrix. For network description and illustration, we calculated the mean and skewness of STTC across all electrode pairs.

In the second part of the package (right side in GUI, compare Supplementary Figures S4C,D), functions are included that are capable of automatically grouping data according to a grouping file (.xlsx) that is provided by the user via the GUI.

Principal component analysis for time series recordings

For PCA, we developed functions in MATLAB that use the well-level data files obtained and saved after running the first part of the routine. Depending on the analysis, these functions created parameter arrays consisting of 19-dimensional parameter vectors (Supplementary Table S2) to quantify cluster separation between the RASopathy model and wild type (wt) using principal components 1 and 2 (PC1/2 analysis). Furthermore, PCA was performed to create vectors of specific parameters describing the features spiking, bursting and network bursting (Supplementary Figure S5B) following up the projection in PC1 with time. After PCA, MATLAB functions grouped well-level data automatically according to genotype or animal/preparation and time. For PC1/2 analysis, PCA score plots visualized genotype clusters in a PC1-PC2 two-dimensional coordinate system. To quantify the separation between the genotype clusters, we followed the suggested quantification method presented in Goodpaster and Kennedy (2011). In brief, we calculated the Mahalanobis distance (Mahalanobis, 1936) resulting in the distance between the centroids between groups considering groups’ data distribution. The calculations were performed by custom-written MATLAB functions that additionally calculated statistical significance with F-statistics.

Electrical evoked activity in stimulation experiments

For an automated analysis of data derived from MEA-based electrical stimulation of neuronal cultures, we developed a software package in MATLAB finally returning the network-wide evoked spiking rate upon STIM normalized to the spontaneous activity (evMFR_norm) (data processing illustrated in Supplementary Figure S6). Based on the initial measurement of spontaneous activity (SA1), preselection of valid wells was performed based on the condition of a minimum number of active electrodes per array (here: 8 electrodes) and with an electrode defined as active whenever its spiking rate exceeded a threshold (herein 0.1 Hz). To calculate evMFR_norm, first, the evoked activity evMFR^’ was determined as the number of spikes within a period of 1 s upon an electrical pulse divided by this time (Supplementary Figure S3B), followed by averaging across all pulses on one electrode during one session. For each electrode, evMFR^’ was then normalized to the spontaneous activity (${\mathrm{M}\mathrm{F}\mathrm{R}}_{\mathrm{S}\mathrm{A}}^{\prime }$) of the same electrode, resulting in ${\mathrm{evMFR}}_{\mathrm{norm}}^{\prime }$. To evaluate the network-wide evoked activity, ${\mathrm{evMFR}}_{\mathrm{norm}}^{\prime }$ was averaged across all array electrodes and across all sessions within STIM, yielding evMFR_norm. The software package returned a data array containing evMFR_norm, as well as the corresponding coefficient of variance and the skewness for each well and STIM. Based on a user-provided Excel file, the package automatically grouped the wells according to genotype.

Data cleaning and selection of valid wells

To generate reliable results and to reduce the variability due to hardly predictable environment and physiological changes, criteria for the selection of valid wells were defined. First, wells were selected according to their activity as measured by the minimum number of active electrodes (if not otherwise stated: 10) determined by the minimum spiking rate (herein >0.1 Hz). For the time series experiments, this selection was performed based on measurements taken at the time of network maturity, which we denote here as the reference time point (DIV 21–24). For disinhibition-induced stimulation experiments, the baseline measurement at the fourth time point was taken as a reference. Wells that did not fulfil the criteria were discarded from the analysis. For time series experiments we implied further quality criteria and discarded wells once they lost their integrity beyond the reference time point. Loss of integrity was defined if both following conditions were given (1) a considerable drop in the number of active electrodes (drop of contributing channels ≥2 channels compared to reference time point) and (2) a clear decrease in neuronal activity (wMFR <70% wMFR_ref). A precise overview of datasets obtained from recordings of spontaneous activity of both genotypes is given in Supplementary Table S5. Detailed description of criteria used for exclusion of data points are listed in Supplementary Table S6. For PCA, wells with missing values over the course of an experiment were discarded.

For stimulation experiments that were more prone to interference due to a measurement duration above 2 h, additional exclusion criteria were applied to ensure that only stable cultures were included in the analyses. First, cultures had to exhibit stability in their spontaneous activity prior to (SA1) and after (SA2) STIM1 (Supplementary Figure S3D). For each well, the MFR and SD of the firing rate were computed, and significant changes between SA1 and SA2 were identified by Welch’s t test to discard wells that were significantly altered. Second of particular importance for the evaluation of the effect of tetanic stimulation on the network evoked activity, we discarded cultures that were pathologically affected by STIM sessions. To evaluate this, we regarded the relation between evMFR_norm upon STIM2 and upon STIM1 and checked for outliers (Supplementary Table S4).

Statistics applied to MEA datasets

For statistical calculations, we used GraphPad Prism (version 8.3.0) Software and MS Excel. The distribution of all data sets was tested by normality and lognormality tests (Supplemetary Table S3 and Supplementary Figures S7, S8). To test significant differences within time series experiments, subsequent PCA and disinhibition-induced stimulation experiments, we performed a mixed-effects analysis with GraphPad Prism. For time series experiments, we assume an increased risk of type two error caused by a systematic effect by time on the whole population potentially superimposing the effect of genotype and therefore performed either an unpaired t test for normally distributed data or Mann–Whitney U tests for nonnormally distributed data without correction for multiple comparisons. The effect size is indicated as r for unpaired student’s t-test (Lakens, 2013) and r for the Mann–Whitney U test according to Wendt (1972) (Figures 1, 2). Standard benchmarks for r fit well with our set of data: <0.3 small effect, 0.3–0.5 moderate effect, and >0.5 strong effect A p-value <0.05 was considered significant.

Figure 1

Figure 1. Spontaneous spiking and bursting activity in neuronal networks from the RASopathy models Ptpn11^D61Y and Kras^V14l. (A,B) Parameters describing neuronal activity as a function of culture age as box plots indicating the median and interquartile range with whiskers extending from 5 to 95% confidence interval. (A) wMFR and interspike interval. (B) Burst duration and wMBR. Parameters are statistically evaluated by a mixed-effects model approach. Below the graphs, p and F values for main and intermediate effects are shown. If the F test reveals significance in genotype, p-values from the Mann–Whitney U test (black) and effective size r according to Wendt (gray, cursive) are shown in the graphs. Stars indicate significant differences between genotypes for a time point: <0.05 (*), <0.01 (**). A data point represents the mean values across the wells related to one animal/preparation (preparation level).

Figure 2

Figure 2. Functional connectivity/synchronicity in networks with RASopathy mutations Ptpn11^D61Y and Kras^V14l. (A) Number of animals and wells analyzed at each culture age cohort for all parameters of spontaneous activity. (B,C) Parameters describing functional connectivity as a function of culture age as box plots indicating the median and interquartile range with whiskers extending from 5 to 95% confidence interval. (B) Values of the number of network bursts. (C) Values of mean STTC and skewness STTC. Parameters are statistically evaluated by a mixed-effects model approach. Below the graph, p and F values for main and intermediate effects are shown. If the F test reveals significance in genotype, in B, p-values from the Mann–Whitney U test (black) and effective size r according to Wendt (gray, cursive) are shown; in C, p-values from the unpaired t-test together with effect size r² are shown in the graphs. Stars indicate significant differences between genotypes for one time point <0.05 (*), <0.01 (**), <0.001 (***), <0.0001 (****). A data point represents the mean values across the wells related to one animal/preparation (preparation level).

The projected data set on PC1 was shown to be approximately normally distributed (Supplementary Figure S9) with equal variances in all groups (Supplementary Figure S10), leading us to perform a mixed-effects model with Sidak’s test to correct for multiple comparisons for increased statistical power. Since concomitant with PCA, normalization of data values centers the data points around zero for each time point, the systematic increase of neuronal activity coming along with network development was filtered out in all groups. Thus, the effect on spiking, bursting and network bursting phenotypes was unveiled, decreasing the probability of error type two. Additionally, for disinhibition-induced stimulation experiments, Sidak’s multiple comparisons test was performed and for effect size Cohen’s d was computed and converted to correlation coefficient r as effect size.

To compare the evoked activity upon STIM, we used the nonparametric Mann–Whitney test due to the nonnormal distribution of the data set for MFR. Additionally, we investigated differences in distribution by the Kolmogorov–Smirnov test. To analyze the effect of tetanic stimulation on networks, we calculated the relation X = $\frac{{\mathrm{evMFR}}_{\text{log}}^{\mathrm{STIM}3}}{{\mathrm{evMFR}}_{\begin{array}{c}\hfill log\hfill \\ \hfill \hfill \end{array}}^{\mathrm{STIM}2}}$ between the logarithmic ${\mathrm{evMFR}}_{\mathrm{norm}}$ after (STIM 3) and prior (STIM 2) tetanic stimulation and checked for alterations between genotypes by Welch’s t test (Welch, 1947). Therefore, we calculated descriptive measures from the independent, individual experiments first, followed by Welch’s t test according to

$t=\frac{{\overline{X}}_{\mathrm{floxed}/\mathrm{W}\mathrm{T}}-{\overline{X}}_{\mathrm{W}\mathrm{T}/\mathrm{W}\mathrm{T}}}{\sqrt{\frac{S_{\mathrm{floxed}/\mathrm{W}\mathrm{T}}^2}{n_{\mathrm{floxed}/\mathrm{W}\mathrm{T}}^2}+\frac{S_{\mathrm{W}\mathrm{T}/\mathrm{W}\mathrm{T}}^2}{n_{\mathrm{W}\mathrm{T}/\mathrm{W}\mathrm{T}}^2}}}$

The mean value $\overline{X}$ is the average of X across all experiments for the indicated genotype. S² is given by the averaged variance of the indicated genotype derived from the individual experiment, and n is the total number of all wells over all experiments for the indicated genotype. To check the effect of tetanus on networks, we additionally performed a one-sample Wilcoxon test on the pooled data set of X from all experiments.

Results

Increase in population-wide neuronal activity in neuronal networks from Ptpn11^D61Ymice

To identify the effect of RASopathy mutations on the functional development of neuronal networks, we performed longitudinal recordings of neuronal network activity in cultured dissociated cortical neurons grown on mwMEA plates. The neurons were derived from the cerebral cortices of newborn mice. As the RASopathy model, we used mice heterozygous for the floxed mutation Ptpn11^D61Y or Kras^V14l as well as for the Emx1-cre allele and as the control, their littermates homozygous for the Ptpn11 or Kras wild-type allele and heterozygous for the Emx1-cre allele. The spontaneous network activity of cortical cultures was recorded multiple times in cultures in a time range of DIV 10 until DIV 33 (time series). For one experiment, cultures were prepared from six to eight newborn mice (P0–P1) littermates. The number of animals and the number of arrays analyzed in each age cohort of the culture for all parameters of spontaneous activity are shown in Figure 2A.

To analyze the parameters describing spontaneous activity in neuronal networks, we averaged the data derived from individual wells across preparation/animal and subsequently grouped the datapoints according to genotype. To assess differences between groups and time points, a mixed-effects model was calculated. Herein, we accepted the violation of the distributional assumption of equal SD within groups due to the robustness of the mixed effects model (Supplementary Figure S8). The global evaluation revealed that all groups reliably exhibited spontaneous activity on DIV 11–13. A gradual increase in population-wide firing activity (assessed as spiking and bursting, Figures 1A,B) and in functional connectivity/network synchronicity (assessed as network bursting and correlation in spiking, Figures 2B,C) was evident between DIV 19 and DIV 25, when all parameters stabilized, presumably marking network maturity. Sample sizes are listed in Figure 2A.

To assess the effect of RASopathy mutations on spiking and bursting behavior, we followed up on four parameters, including wMFR, interspike interval, burst duration and wMBR (Figure 1). Ptpn11^D61Y networks were significantly different in most parameters compared to their control according to the F test in the mixed effects model. The spiking activity was increased as reflected by significantly higher values in wMFR at the early time points (DIV 11–13 and DIV 15–16) and in a decreased interspike interval reaching significance throughout all observation time points except on DIV 30–33 (Figure 1A, left column). Furthermore, Ptpn11^D61Y networks showed increased bursting behavior compared to their control (Figure 1B, left column). The appearance of spontaneous bursts during development is important for the formation of neuronal circuits (Kirkby et al., 2013). Ptpn11^D61Y networks exhibited significantly higher wMBR on DIV 11–13, DIV 15–16 and DIV 19–20. However, the burst duration remained unchanged. In contrast, except the significantly increased wMBR at DIV 17–19 we detected no significant differences in any parameter describing spiking and bursting behavior in Kras^V14l. However, visualization of the data indicates a conspicuous trend of higher median values in parameters reflecting neuronal activity (wMFR, wMBR, burst duration) (Figures 1A,B, right column).

A fundamental feature of neuronal networks is their functional connectivity/network synchronicity, which describes the impact of the activity of other neurons within the same network on the activity of a given neuron (Makarov et al., 2005). To assess functional connectivity during network development, we calculated the number of network bursts (Figure 2B) and determined the correlation between all pairs of electrodes during NBs measured as STTC (Figure 2C). Ptpn11^D61Y networks globally demonstrated significantly higher functional connectivity than their control groups, as revealed by a mixed effects model in all observed parameters (Figures 2B,C, left side). In Ptpn11^D61Y networks, the number of NBs was increased at early time points (DIV 11–13, DIV 15/16 and DIV 19/20), but this effect was abolished at later time points. Moreover, compared to controls, Ptpn11^D61Y networks indicated significant differences in the mean and skewness of STTC throughout the whole observation time, as tested by multiple comparison tests (Figure 2C, left side). In contrast, compared to the control group, Kras^V14l showed no alterations in the number of NBs, in the mean STTC and in the STTC skewness (Figures 2B,C, right side). However, a significant global difference was calculated for the skewness of STTC according to the F test on the main effect genotype on DIV 15–16. In general, we identified elevated predicted median values in the Kras^V14l networks in most parameters; however, the differences were not significant, likely due to large variability, which is typical of MEA-derived data.

Principal component analysis of MEA parameters reveals differences in genotype clusters and identifies differences in spiking, bursting and network bursting in both RASopathy models

Analysis of most individual parameters from the longitudinal recordings in Ptpn11^D61Y neurons revealed higher effect sizes in younger cultures in several parameters (wMFR, ISI, wMBR, no. of network bursts) that were reduced at later time points. Similar tendency (i.e., large effect size r) was observed in Kras^V14l neurons, although it did not reach significance. To increase the interpretability of data sets and the power of statistical testing and to allow simpler visualization, we introduced PCA, which reduces data dimensionality. We performed PCA including 19 parameters obtained from a custom-written MATLAB program (Supplementary Table S2). We conducted the PCA at the well level. In contrast to the preparation level, we yielded a higher number of data points, resulting in a lower probability of error type two but decreased power.

We analyzed all age cohorts from the longitudinal recordings and visualized the transformed data points into a PC1 and PC2 coordination system (Figures 3A, 4A). Sample numbers are shown in Figures 3D, 4D. In both RASopathy disease models, the respective clusters and their controls appear to converge with culture age, also indicated by the visualized distance between the cluster centroids in the graphs. To quantify the separation between genotype clusters, we calculated the Mahalanobis distance that accounts for the correlation of the variables (Mahalanobis, 1936). Statistical significance was evaluated by Hotelling’s two-sample T² related to an F value as introduced earlier (Goodpaster and Kennedy, 2011). For both RASopathy model networks, the Mahalanobis distance to the control group turned out to be strongly significant at the early stages but dramatically decreased with culture age (Figures 3B, 4B). At the end of the observation period, Ptpn11^D61Y networks still differed significantly from their control group according to the F-statistic on DIV 30–33, while for Kras^V14l, the clusters converged more closely and did not differ significantly from DIV 24/25 on. The reliability of the results is underpinned by the cumulative fraction of variance of PC1 and PC2 explaining approximately 70% of the variation in both data sets (Figures 3C, 4C). This approach confirmed a developmental phenotype in spontaneous network activity in both RASopathy disease models.

Figure 3

Figure 3. PCA of the Ptpn11^D61Y 19-dimensional parameter vector. (A) Principal component PC1 versus PC2 scores plotted as a function of culture age (DIV). A solid line is drawn between the centroids of the control cluster (triangle) and mutant cluster (rectangle). Each color-coded data point represents a recording from one well. (B) Cluster separation in dependence of culture age was measured by Mahalanobis distance (green dotted line), and the F value shows statistical significance (black solid line) (F_critical = 3.1). (C) Fraction of variance explained by each PC summed over all PCs. The cumulative fraction of variance is indicated for each evaluated culture age. (D) Sample size is given per culture age as the number of animals and wells.

Figure 4

Figure 4. PCA of the Kras^V14l 19-dimensional parameter vector. (A) Principal component PC1 versus PC2 scores plotted as a function of culture age (DIV). A solid line is drawn between the centroids of the control cluster (triangle) and mutant cluster (rectangle). Each color-coded data point represents a recording from one well. (B) Cluster separation in dependence of culture age was measured by Mahalanobis distance (green dotted line) between the clusters of genotypes, and the F value shows statistical significance (black solid line) (F_critical = 3.1). (C) Fraction of variance explained by each PC summed over all PCs. The cumulative fraction of variance is indicated for each evaluated culture age. (D) Sample size is given per culture age as the number of animals and wells.

To further identify at which level spontaneous network activity is particularly affected in RASopathy models, we performed a further PCA that specifically combines parameters describing the features of spiking, bursting and network bursting behavior (Supplementary Figure S5B). Therefore, we formed parameter vectors from the extracted parameters from the spontaneous activity analysis. The values were logarithmized for nonnormally distributed data. This resulted in a 4-dimensional vector for spiking (log₁₀(wMFR), log₁₀(wtMFR), log₁₀(interspike interval), and number of contributing electrodes), a 3-dimensional vector for bursting (log₁₀(wMBR), log₁₀(burst duration), and log₁₀(number of spikes per burst)) and a 4-dimensional parameter vector for network bursting (number of network bursts, mean NB duration, number of spikes per NB, and number of contributing channels). The PCA resulted in a component projection on PC1 for each age cohort (Figure 5). We visualized the projected parameter vector on PC1 for each neuronal age cohort (Figures 5A–C). We were able to identify a changed phenotype in spiking, bursting and network bursting for both RASopathy models compared to wt.

Figure 5

Figure 5. PCA of parameter vectors related to spiking, bursting and network bursting from Ptpn11^D61Y and Kras^V14l. (A–C) PC1 projection as a function of culture age as box plots indicating the median and interquartile range with whiskers extending from 5 to 95% confidence interval. (A) Quantification of spiking behavior calculated as a projection of a 4-dimensional parameter vector (number of spiking electrodes, wMFR, wtMFR, interspike interval). (B) Quantification of bursting behavior calculated as the projection of a 3-dimensional parameter vector (wMBR, bursting duration, no. of spikes per burst). (C) Quantification of network bursting calculated as a projection of a 4-dimensional parameter vector [number of network bursts (NB), mean NB duration, number of spikes per NB, mean number of contributing channels per NB]. Statistical significance was assessed by a mixed-effects model approach. p-values (in black) for each comparison were obtained by Sidak’s multiple comparisons test, and effect size (r², in gray, cursive) was derived from the unpaired t-test. Stars indicate significant differences in median values between genotypes for time points: 0.05 (*), 0.01 (**), 0.001 (***), <0.0001 (****).

On the spiking level, the mixed effects model revealed a strongly significant effect of genotype for both RASoapthy models (Figure 5A). Multiple comparisons testing showed the strong differences primarily at early time points (Ptpn11^D61Y: DIV 11–13, DIV 15/16, DIV 19/20; Kras^V14l: DIV 15/16, DIV 17–19, DIV 20–22). Similarly, for network bursting, the mixed effects model indicated a significant effect of genotype for both mutations, with differences on DIV 11–13 and on DIV 19/20 in Ptpn11^D61Y and on DIV 15/16 and on DIV 17–19 in Kras^V14l, respectively (Figure 5C). In contrast, as assessed by multiple comparisons testing, for bursting, Ptpn11^D61Y showed significant differences at latest time point (DIV 30–33), when networks are expected to be mature, while Kras^V14l reveals significant differences at DIV 17–19 and DIV24/25. Interestingly, the mixed effects model identified a strong interaction effect over time and between genotypes in all parameter vectors describing the network activity for both RASopathy models (genotype × time, indicated in all graphs in Figure 5) indicating differences in the functional maturation of neuronal networks.

Overall, we observed strong significant differences in spiking for both RASopathy models, that were more pronounced during early development and shades out during network maturation. The changes in bursting and network bursting were also significant, but observed only at specific timepoints.

Pharmacological disinhibition has a specific effect on network bursting in Ptpn11^D61Y

Our experiments revealed a change in neuronal activity and functional connectivity during network development in both RASopathy models. Interestingly, the effect of mutations on network characteristics is attenuated over time, pointing to possible adaptation. To assess the possible contribution of a shift in excitatory/inhibitory (E/I) balance, as previously suggested in Papale et al. (2017), we decided to test the effect of pharmacological interference with inhibitory tone on spontaneous network activity. Therefore, we applied bic, a potent antagonist of the GABA_A receptor, to cultured cortical neurons (DIV 33) from Ptpn11^{D61Y mice} and their wild-type littermates (Figure 6A). As expected, the bic treatment resulted in a reliable increase in spiking activity reflected in increased wMFR (Figures 6B,C), which was statistically confirmed by one sample Wilcoxon test for each well (Figure 6C). While spiking and bursting behavior did not differ between genotypes upon bic treatment (Figure 6C; Supplementary Figure S11), we detected elevated network bursting in disinhibited Ptpn11^D61Y networks compared to the control. The mixed effects model revealed a global interaction effect of time and genotype (genotype x time p < 0.0001), leading us to perform Sidak’s multicomparison tests resulting in significant differences from minute 50 post bic (Figures 6C,D). For an overall conclusion of the main effect genotype, the mixed effects model returns a p-value of.06 that is near significance. This result revealed changes in the inhibitory system in Ptpn11^D61Y networks, which likely developed to counteract increased spontaneous network activity.

Figure 6

Figure 6. Effect of disinhibition on neuronal network activity in Ptpn11^D61Y on DIV 33. (A) Number of animals (from 3 independent experiments) and wells used for this analysis. (B) Representative raster plots showing spatiotemporal spiking activity of control and Ptpn11^D61Y networks prior to and 30 min after bicuculline application (bic). (C,D) Line graphs of wMFR (C) and MNBR (D) as a function of time prior to and after application of bic showing an elevated spiking upon treatment in both genotypes. Values are represented in percent relative to a baseline generated from values recorded over a time interval of 20 min prior to bic treatment, averaged and set to 100%. A significant increase against baseline was statistically tested by one-sample Wilcoxon. In (D), the interaction effect between time and genotype was calculated by a mixed effects model approach. Sidak’s multiple comparison tests were used to determine significance between genotypes. Stars indicate significant differences <0.05 (*). To assess the possible contribution of a shift in excitatory/inhibitory (E/I) balance, as previously suggested in Papale et al. (2017), we decided to test the effect of pharmacological interference with inhibitory tone on.

Fewer GABAergic synapses but higher expression of vesicular transporters for GABA and glutamate were detected in the networks with the Ptpn11^D61Y mutation

To examine the altered inhibitory tone on a synaptic level, we visualized and quantified glutamatergic and GABAergic boutons in Ptpn11^D61Y networks and their control. Therefore, neuronal cultures were immunolabeled with specific antibodies against vesicular GABA transporter (VGAT) or vesicular glutamate transporter (VGLUT1), together with antibodies against synaptophysin (Syp), to visualize all synapses (Figures 7A,B). We used our in-house developed MATLAB-based program SynEval to obtain Syp-positive puncta (ROIs) and to sort them as VGAT- or VGLUT-positive in an automatized batch mode, which reduces human bias and time effort. The fraction of glutamatergic boutons was significantly higher in the Ptpn11^D61Y group than in the control group (Figure 7D), and accordingly, the fraction of GABAergic boutons decreased (Figure 7C). Interestingly, we noticed an increased synaptic IFI for VGAT or VGLUT, indicating increased synaptic abundance of these neurotransmitter transporters (Figures 7C,D). These findings further point toward changes in the excitation-inhibition balance in Ptpn11^D61Y networks.

Figure 7

Figure 7. Immunolabeling of excitatory and inhibitory synapses in neuronal networks from Ptpn11^D61Y mice. (A) Representative images of neurons (DIV 21) labeled using antibodies against VGAT (green) and Syp (red). (B) Representative images of neurons labeled using antibodies against VGLUT1 (green) and synaptophysin (Syp, red). (C) Quantification of VGAT IFI in Syp-positive synaptic puncta and of the fraction of VGAT-positive synapses. (D) Quantification of VGLUT1 IFI in Syp-positive synaptic puncta and of the fraction of VGLUT1-positive synapses. (C,D) Sample size is given as the number of analysed images. Statistical significance was determined by the Mann–Whitney test. Data points correspond to individual analyzed visual fields derived from two independent experiments. The scale bar is 5 μM.

Altered evoked activity in Ptpn11^D61Y cells probed by electrical stimulation

The longitudinal analysis indicated a normalization of spontaneous network activity during maturation in the neuronal networks expressing RASopathy mutations. To assess whether there are functional defects in the processing of information in these apparently “normal” networks, we analyzed network activity evoked by the delivery of electrical stimulation. We recorded network activity evoked by single-pulse stimulation (test stimulus, STIM) delivered from different electrodes in Ptpn11^D61Y networks grown for 33 days (DIV 33) (Supplementary Figures S3A,D). Of note, we solely analyzed robust and stable networks, which were not affected by the application of STIM itself, i.e., the activity evoked by STIM1 was the same as the activity evoked by STIM2 (cf. Material and Methods, Data cleaning and selection of valid wells). The data points plotted in Figure 8A display the evMFR_norm representing the evoked MFR induced by STIM at the electrode level and normalized to recorded spontaneous activity on the same electrode, averaged across an array (for detailed information on calculation, see Supplementary Figure S6). Interestingly, evMFR_norm was found to be significantly decreased in Ptpn11^D61Y networks compared to wt (Figure 8A, upper plot). The graph shows the data for STIM2, and similar results were obtained for STIM1 (Supplementary Figure S12). The Kolmogorov–Smirnov test confirmed a different distribution of data between both genotypes.

Figure 8

Figure 8. Assessment of evoked network activity in Ptpn11^D61Y neurons. (A) Evoked mean firing rate (evMFR_normed) before and after TET stimulation. Black lines indicate median values, and data points correspond to one well. Significance was tested by the Mann–Whitney (MW) test to check differences in mean values and the Kolmogorov–Smirnov (KS) test to check differences in data distributions. (B) Relative change (ratio) of evMFR_normed recorded before and after application of TET are shown from three independent experiments (Exp1–Exp3). Numbers indicate the group sizes for the experiment (animals/wells). Whiskers indicate SEM. The p-value was calculated by Welch’s t test; the number of wells included for the corresponding genotype is color-coded for each experiment. (C) Analysis of data pooled from all 3 experiments. Black lines indicate median values, and one data point refers to one well. The one-sample Wilcoxon test (WT) was tested against a hypothetical value of 1 (indicating no change in evMFR_normed caused by TET), and the Mann–Whitney test (MW) was used to calculate significant differences between genotypes.

Effect of high-frequency electrical stimulation on evoked activity in Ptpn11^D61Y mice

We observed that evoked neuronal activity upon STIM is affected by Ptpn11^D61Y mutations and that excitability is dampened. Next, to understand whether neuronal plasticity might also be affected in Ptpn11^D61Y networks, we tested the effect of tetanic stimulus (TET) on evoked synaptic activity. The stimulation protocol described in Supplementary Figure S3C was previously suggested to be suitable for analysis of neuronal plasticity since it affected evoked activity in networks (Chiappalone et al., 2008). The TET stimulus applied on Ptpn11^D61Y networks induced different effects compared to controls in all three independent experiments individually plotted in Figure 8B. Calculating the relative change (ratio) between post and pre TET logarithmic (log₁₀)evMFR_norm showed that the mean value in Ptpn11^D61Y was higher than that in the control group for each experiment (Exp1-Exp3), with an overall significant difference confirmed by Welch’s t test (Mann–Whitney test for individual experiments: Exp1: p = 0.03, Exp2: ns, Exp3: ns). More precisely, the findings reveal that TET resulted in a significant increase in evMFR_norm in Ptpn11^D61Y networks while not inducing an effect in controls (Figure 8C). That was shown by testing the pooled data set against the hypothetical value 1 (no change between pre and post TET logarithmic ${\mathrm{evMFR}}_{\mathrm{norm}}$ ) by a one-sample Wilcoxon test. Interestingly, the difference in evMFR_norm, demonstrated between Ptpn11^D61Y and control pre TET, was diminished post TET (Figure 8A, bottom plot). According to the Kolmogorov–Smirnov test, the difference in distribution also failed significance after TET. Indeed, in contrast to the control, Ptpn11^D61Y demonstrated a shift between STIM2 and STIM3 curves in their relative frequency plots, particularly for higher ${\mathrm{evMFR}}_{\mathrm{norm}}$ , pointing to an altered distribution in the data set (Supplementary Figures S13B,C).

Discussion

RASopathies are a group of rare genetic disorders with various symptoms, including heart defects, skeletal anomalies, increased risk for tumors and variable degrees of neurocognitive deficits. Clinically, the severity of neurocognitive impairments in RASopathies often varies over the lifetime and often becomes milder with the age of the patients. Neurodevelopmental phenotypes were also described in RASopathy mouse models, where the involvement of maladaptive changes in cortical networks was suggested Papale et al. (2017).

To investigate whether similar maladaptation occurs in neuronal networks on a chip, we cultured cortical neurons carrying two distinct RASopathy mutations on mwMEAs and recorded their population-wide spontaneous network activity to monitor network establishment and maturation during a time period of 5 weeks. To facilitate the evaluation of large data sets acquired longitudinally with a high sampling rate (12.5 kHz), we developed and implemented an analysis pipeline programmed in MATLAB. From the list of spikes per electrode, provided by acquisition software, it computes a custom set of parameters describing various features of neuronal spontaneous network activity, including spiking, bursting, network bursting and correlation in spiking. This custom-written software package permits fast and unbiased analysis of large MEA data sets while being flexibly applicable for a wide set of experimental designs in general. The provision of a GUI renders the package user-friendly. Furthermore, we developed an add-on package of MATLAB-based functions that implements PCA to reduce the complexity of analysis outcomes by projecting the resulting parameters into a two-dimensional system. The package allows the calculation of Mahalanobis distance and statistical testing to mathematically compare the network behavior between two groups derived in our case from mutant and wild-type animals and determine if these constitute significantly distinguishable populations. In the next step, it also computes whether the difference between network clusters can contribute to one specific network function feature (spiking, bursting and network synchronicity) or whether changes in multiple features contribute to the clusters’ differences. This approach turned out to be suitable to uncover effects in the MEA data that are inherently characterized by high complexity and variability.

On the basis of MEA-recorded data and supported by the MATLAB analysis pipeline, we identified a convergent developmental phenotype in both RASopathy models. Multiple parameters describing neuronal activity and functional connectivity were significantly different in neurons with the Ptpn11^D61Y mutation. The differences were more pronounced in younger networks and faded out as networks matured. In contrast, the neurons with the Kras^V14l mutation showed only a tendency toward similar changes in spontaneous network activity parameters as seen in Ptpn11^D61Y neurons. Being aware of the high inherent variance of MEA data, we performed PCA to obtain an overall description of network parameters for both models. This analysis revealed a significant separation between clusters representing neurons with the Ptpn11^D61Y or Kras^V14l mutation and their respective controls. Moreover, for both mutations, the cluster separation from control networks decreased throughout development, even though for Ptpn11^D61Y the difference in the individual parameters mean STTC and skewness of STTC remains significant at DIV 30. Consecutive analysis identified a major contribution of changes in spiking and network bursting to the PCA cluster separation. Thus, based on the analysis of spontaneous network activity according to PCA, we have clearly demonstrated developmental phenotypes for both mutations, confirming that the effect of RASopathy mutations on neuronal network function converges. It is interesting to note that the strength of mutations’ effect on network function correlates with the strength of their activation assessed biochemically in previous studies. While the Noonan syndrome associated Kras^V14l exhibits attenuated biochemical activation compared to common cancer associated KRAS alleles (Schubbert et al., 2007), strong activation was shown for Ptpn11^D61Y, which is frequently linked to juvenile myelomonocytic leukemia (JMML) in patients (Chan et al., 2009). Still, both mouse models used in this study are appropriate to investigate the neurodevelopmental effect of RASopathies. The mice with cortex specific expression of Ptpn11^D61Y display lower exploratory activity, reduced memory specificity, and alteration of neuronal activity-induced gene expression (Altmuller et al., 2017). Similarly, the mice with neuron-specific expression of oncogenic KrasG12V reliably model Noonan-syndrome linked neurocognitive symptoms (Papale et al., 2017). Of note, the cortical cultures used in this study contain next to neurons also astrocytes and to lesser extent also oligodendrocytes. Non-neuronal cells affect the functional maturation and network activity of neuronal networks on chip as well as in vivo (Enright et al., 2020). It is therefore feasible that expression of RASopathy mutations in astrocytes as it occurs in our model system as well as in patients contribute to the phenotypes observed in this study.

The normalization of network activity parameters in mature networks suggests that compensatory mechanisms develop to counteract the increased network activity observed in juvenile networks. A shift in the E/I balance might represent one such mechanism. Indeed, the balance between E/I within neuronal circuits is crucial for correct brain development and function, and its disturbances contribute to numerous neurodevelopmental disorders, such as autism spectrum disorders and epilepsy, which are common comorbidities in more severe RASopathies, including CFC or CS (Adviento et al., 2014; Pierpont et al., 2022). Moreover, the recent characterization of a RASopathy animal model expressing the strongly activating allele Kras^G12V indicated that increased inhibitory drive contributes to the electrophysiological and behavioral abnormalities in this model Papale et al. (2017). In line with these assumptions, we demonstrated increased network bursting in mature Ptpn11^D61Y networks upon pharmacological disinhibition of GABA receptors, which indicates an increase in the inhibitory tone. However, the immunocytochemical analysis revealed an increased proportion of excitatory and decreased proportion of inhibitory synapses in Ptpn11^D61Y networks, but an increased expression of vesicular transporters for both glutamate and GABA. Thus, although it is difficult to directly correlate morphology and function, our morphological findings confirm that specific network and synaptic adaptation occur in Ptpn11^D61Y neuronal networks. Further, we found that the application of electrical test stimuli elicited a reduced evoked response in the Ptpn11^D61Y networks. The attenuated response to stimulation points toward maladaptive dampening of network activity and is well compatible with a compensatory increase in inhibitory tone. An application of TET stimulation that had no effect on evoked responses in the control had a significant impact on responses in neurons with the Ptpn11^D61Y mutation. Interestingly, TET stimulation statistically restored the affected evoked activity in the RASopathy model. However, the change in the distribution of data values between pre- and posttetanus stimulations indicates that there is much higher variability in the responses of individual tested networks expressing Ptpn11^D61Y mutations compared to controls. Taking together, while several results of our study point to a compensatory increase in inhibitory tone, the exact mechanisms involved in RASopathy-induced network adaptation require further investigation. A better understanding of the underlying cellular and molecular mechanisms will facilitate therapeutic approaches to disrupt the maladaptive circuits and develop targeted treatments for RASopathies and their comorbidities.

In summary, our study reveals that Ptpn11^D61Y and Kras^V14l mutations have convergent effects on neuronal network activity. We observed more severe differences in network characteristics in juvenile networks that were compensated in the late stages of network maturation, which is in good agreement with clinical data. Despite the apparent normalization of spontaneous network activity, we uncovered significant differences in evoked responses in mature networks, which are likely due to changes in E/I balance. Finally, the cell on-a-chip approach developed here is well scalable, accessible for treatments and manipulations and adaptable for stem cell technology-derived humanized disease modeling. Therefore, we assume that the analysis workflow developed here could be implemented for drug screening and testing platforms also for unrelated brain diseases.

Data availability statement

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.

Ethics statement

The animal study was approved by Sachgebiet Tierschutz, FAU Erlangen. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

E-MW: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. DG: Investigation, Writing – review & editing. AP: Investigation, Visualization, Writing – review & editing. VH: Investigation, Writing – review & editing. MZ: Funding acquisition, Resources, Writing – review & editing. AF: Conceptualization, Funding acquisition, Methodology, Resources, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The study was funded by BMBF Project GeNeRARe: German Research Network for RASopathies (01GM1902B). This work was cofounded by EURAS (Grant Agreement No. 101080580) and by the ELAN Program (Grant Number P112) of the medical faculty of the Friedrich-Alexander-Universität Erlangen-Nürnberg.

Acknowledgments

We are deeply grateful to all those who contributed to the success of this research project. First, we would like to thank our cooperation partner from the stem cell biology department at the University Hospital Erlangen, particularly Prof. Dr. Beate Winner, Dr. Tania Rizo and Michaela Farrell for great technical support. Furthermore, we thank the group of Prof. Dr. Thomas Winkler at the Department of Biology at Friedrich-Alexander University of Erlangen for their support.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2024.1388409/full#supplementary-material

References

Adviento, B., Corbin, I. L., Widjaja, F., Desachy, G., Enrique, N., Rosser, T., et al. (2014). Autism traits in the RASopathies. J. Med. Genet. 51, 10–20. doi: 10.1136/jmedgenet-2013-101951

PubMed Abstract | Crossref Full Text | Google Scholar

Alfieri, P., Piccini, G., Caciolo, C., Perrino, F., Gambardella, M. L., Mallardi, M., et al. (2014). Behavioral profile in RASopathies. Am. J. Med. Genet. A 164, 934–942. doi: 10.1002/ajmg.a.36374

PubMed Abstract | Crossref Full Text | Google Scholar

Altmuller, F., Pothula, S., Annamneedi, A., Nakhaei-Rad, S., Montenegro-Venegas, C., Pina-Fernandez, E., et al. (2017). Aberrant neuronal activity-induced signaling and gene expression in a mouse model of RASopathy. PLoS Genet, 13, e1006684.

Google Scholar

Altmuller, F., Pothula, S., Annamneedi, A., Nakhaei-Rad, S., Montenegro-Venegas, C., Pina-Fernández, E., et al. (2017). Aberrant neuronal activity-induced signaling and gene expression in a mouse model of RASopathy. PLoS Genet. 13:e1006684. doi: 10.1371/journal.pgen.1006684

PubMed Abstract | Crossref Full Text | Google Scholar

Anni, D., Weiss, E. M., Guhathakurta, D., Akdas, Y. E., Klueva, J., Zeitler, S., et al. (2021). Abeta1-16 controls synaptic vesicle pools at excitatory synapses via cholinergic modulation of synapsin phosphorylation. Cell. Mol. Life Sci. 78, 4973–4992. doi: 10.1007/s00018-021-03835-5

Crossref Full Text | Google Scholar

Araki, T., Chan, G., Newbigging, S., Morikawa, L., Bronson, R. T., and Neel, B. G. (2009). Noonan syndrome cardiac defects are caused by PTPN11 acting in endocardium to enhance endocardial-mesenchymal transformation. Proc Natl Acad Sci U S A, 106, 4736–41.

Google Scholar

Armour, C. M., and Allanson, J. E. (2008). Further delineation of cardio-facio-cutaneous syndrome: clinical features of 38 individuals with proven mutations. J. Med. Genet. 45, 249–254. doi: 10.1136/jmg.2007.054460

PubMed Abstract | Crossref Full Text | Google Scholar

Bang, S., Hwang, K. S., Jeong, S., Cho, I. J., Choi, N., Kim, J., et al. (2021). Engineered neural circuits for modeling brain physiology and neuropathology. Acta Biomater. 132, 379–400. doi: 10.1016/j.actbio.2021.06.024

Crossref Full Text | Google Scholar

Cesarini, L., Alfieri, P., Pantaleoni, F., Vasta, I., Cerutti, M., Petrangeli, V., et al. (2009). Cognitive profile of disorders associated with dysregulation of the RAS/MAPK signaling cascade. Am. J. Med. Genet. A 149A, 140–146. doi: 10.1002/ajmg.a.32488

Crossref Full Text | Google Scholar

Chan, G., Kalaitzidis, D., Usenko, T., Kutok, J. L., Yang, W., Mohi, M. G., et al. (2009). Leukemogenic Ptpn11 causes fatal myeloproliferative disorder via cell-autonomous effects on multiple stages of hematopoiesis. Blood 113, 4414–4424. doi: 10.1182/blood-2008-10-182626

PubMed Abstract | Crossref Full Text | Google Scholar

Chiappalone, M., Massobrio, P., and Martinoia, S. (2008). Network plasticity in cortical assemblies. Eur. J. Neurosci. 28, 221–237. doi: 10.1111/j.1460-9568.2008.06259.x

Crossref Full Text | Google Scholar

Costa, R. M., Federov, N. B., Kogan, J. H., Murphy, G. G., Stern, J., Ohno, M., et al. (2002). Mechanism for the learning deficits in a mouse model of neurofibromatosis type 1. Nature 415, 526–530. doi: 10.1038/nature711

PubMed Abstract | Crossref Full Text | Google Scholar

Cui, Y., Costa, R. M., Murphy, G. G., Elgersma, Y., Zhu, Y., Gutmann, D. H., et al. (2008). Neurofibromin regulation of ERK signaling modulates GABA release and learning. Cell 135, 549–560. doi: 10.1016/j.cell.2008.09.060

Crossref Full Text | Google Scholar

Cutts, C. S., and Eglen, S. J. (2014). Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves. J. Neurosci. 34, 14288–14303. doi: 10.1523/JNEUROSCI.2767-14.2014

Crossref Full Text | Google Scholar

Enright, H. A., Lam, D., Sebastian, A., Sales, A. P., Cadena, J., Hum, N. R., et al. (2020). Functional and transcriptional characterization of complex neuronal co-cultures. Sci. Rep. 10:11007. doi: 10.1038/s41598-020-67691-2

PubMed Abstract | Crossref Full Text | Google Scholar

Goodpaster, A. M., and Kennedy, M. A. (2011). Quantification and statistical significance analysis of group separation in NMR-based metabonomics studies. Chemometr. Intell. Lab. Syst. 109, 162–170. doi: 10.1016/j.chemolab.2011.08.009

PubMed Abstract | Crossref Full Text | Google Scholar

Gremer, L., Merbitz-Zahradnik, T., Dvorsky, R., Cirstea, I. C., Kratz, C. P., Zenker, M., et al. (2011). Germline KRAS mutations cause aberrant biochemical and physical properties leading to developmental disorders. Hum. Mutat. 32, 33–43. doi: 10.1002/humu.21377

PubMed Abstract | Crossref Full Text | Google Scholar

Guhathakurta, D., Akdaş, E. Y., Fejtová, A., and Weiss, E. M. (2022). Development and application of automatized routines for optical analysis of synaptic activity evoked by chemical and electrical stimulation. Front. Bioinform. 2:814081. doi: 10.3389/fbinf.2022.814081

PubMed Abstract | Crossref Full Text | Google Scholar

Hernandez-Porras, I., Fabbiano, S., Schuhmacher, A. J., Aicher, A., Cañamero, M., Cámara, J. A., et al. (2014). K-RasV14I recapitulates Noonan syndrome in mice. Proc. Natl. Acad. Sci. U.S.A. 111, 16395–16400. doi: 10.1073/pnas.1418126111

PubMed Abstract | Crossref Full Text | Google Scholar

Kang, M., and Lee, Y. S. (2019). The impact of RASopathy-associated mutations on CNS development in mice and humans. Mol. Brain 12:96. doi: 10.1186/s13041-019-0517-5

PubMed Abstract | Crossref Full Text | Google Scholar

Kim, Y. E., and Baek, S. T. (2019). Neurodevelopmental aspects of RASopathies. Mol. Cells 42, 441–447. doi: 10.14348/molcells.2019.0037

Crossref Full Text | Google Scholar

Kirkby, L. A., Sack, G. S., Firl, A., and Feller, M. B. (2013). A role for correlated spontaneous activity in the assembly of neural circuits. Neuron 80, 1129–1144. doi: 10.1016/j.neuron.2013.10.030

PubMed Abstract | Crossref Full Text | Google Scholar

Krab, L. C., Aarsen, F. K., de Goede-Bolder, A., Catsman-Berrevoets, C. E., Arts, W. F., Moll, H. A., et al. (2008). Impact of neurofibromatosis type 1 on school performance. J. Child Neurol. 23, 1002–1010. doi: 10.1177/0883073808316366

PubMed Abstract | Crossref Full Text | Google Scholar

Lakens, D. (2013). Calculating and reporting effect sizes to facilitate cumulative science: a practical primer for t-tests and ANOVAs. Front. Psychol. 4:863. doi: 10.3389/fpsyg.2013.00863

Crossref Full Text | Google Scholar

Lee, Y. S., Ehninger, D., Zhou, M., Oh, J. Y., Kang, M., Kwak, C., et al. (2014). Mechanism and treatment for learning and memory deficits in mouse models of Noonan syndrome. Nat. Neurosci. 17, 1736–1743. doi: 10.1038/nn.3863

PubMed Abstract | Crossref Full Text | Google Scholar

Mahalanobis, P. C. (1936). On the generalised distance in statistics. Sankhya A 80, 1–7. doi: 10.1007/s13171-019-00164-5

Crossref Full Text | Google Scholar

Makarov, V. A., Panetsos, F., and de Feo, O. (2005). A method for determining neural connectivity and inferring the underlying network dynamics using extracellular spike recordings. J. Neurosci. Methods 144, 265–279. doi: 10.1016/j.jneumeth.2004.11.013

Crossref Full Text | Google Scholar

Mautner, V. F., Kluwe, L., Thakker, S. D., and Leark, R. A. (2002). Treatment of ADHD in neurofibromatosis type 1. Dev. Med. Child Neurol. 44, 164–170. doi: 10.1111/j.1469-8749.2002.tb00780.x

Crossref Full Text | Google Scholar

Molina, J. R., and Adjei, A. A. (2006). The Ras/Raf/MAPK pathway. J. Thorac. Oncol. 1, 7–9. doi: 10.1016/S1556-0864(15)31506-9

Crossref Full Text | Google Scholar

Papale, A., Disa, R., Menna, E., Cerovic, M., Solari, N., Hardingham, N., et al. (2017). Severe Intellectual Disability and Enhanced Gamma-Aminobutyric Acidergic Synaptogenesis in a Novel Model of Rare RASopathies. Biol Psychiatry, 81, 179–192. doi: 10.1016/j.biopsych.2016.06.016

Crossref Full Text | Google Scholar

Payne, J. M., Hyman, S. L., Shores, E. A., and North, K. N. (2011). Assessment of executive function and attention in children with neurofibromatosis type 1: relationships between cognitive measures and real-world behavior. Child Neuropsychol. 17, 313–329. doi: 10.1080/09297049.2010.542746

Crossref Full Text | Google Scholar

Pierpont, E. I., Kenney-Jung, D. L., Shanley, R., Zatkalik, A. L., Whitmarsh, A. E., Kroening, S. J., et al. (2022). Neurologic and neurodevelopmental complications in cardiofaciocutaneous syndrome are associated with genotype: a multinational cohort study. Genet. Med. 24, 1556–1566. doi: 10.1016/j.gim.2022.04.004

PubMed Abstract | Crossref Full Text | Google Scholar

Pierpont, E. I., Tworog-Dube, E., and Roberts, A. E. (2015). Attention skills and executive functioning in children with Noonan syndrome and their unaffected siblings. Dev. Med. Child Neurol. 57, 385–392. doi: 10.1111/dmcn.12621

PubMed Abstract | Crossref Full Text | Google Scholar

Roberts, A., Allanson, J., Jadico, S. K., Kavamura, M. I., Noonan, J., Opitz, J. M., et al. (2006). The cardiofaciocutaneous syndrome. J. Med. Genet. 43, 833–842. doi: 10.1136/jmg.2006.042796

PubMed Abstract | Crossref Full Text | Google Scholar

Roelofs, R. L., Janssen, N., Wingbermühle, E., Kessels, R. P. C., and Egger, J. I. M. (2016). Intellectual development in Noonan syndrome: a longitudinal study. Brain Behav. 6:e00479. doi: 10.1002/brb3.479

PubMed Abstract | Crossref Full Text | Google Scholar

Samuels, I. S., Saitta, S. C., and Landreth, G. E. (2009). MAP’ing CNS development and cognition: an ERKsome process. Neuron 61, 160–167. doi: 10.1016/j.neuron.2009.01.001

PubMed Abstract | Crossref Full Text | Google Scholar

Schubbert, S., Bollag, G., Lyubynska, N., Nguyen, H., Kratz, C. P., Zenker, M., et al. (2007). Biochemical and functional characterization of germ line KRAS mutations. Mol Cell Biol, 27, 7765–70.

Google Scholar

Schubbert, S., Zenker, M., Rowe, S. L., Böll, S., Klein, C., Bollag, G., et al. (2006). Germline KRAS mutations cause Noonan syndrome. Nat. Genet. 38, 331–336. doi: 10.1038/ng1748

Crossref Full Text | Google Scholar

Schulz, A. L., Albrecht, B., Arici, C., van der Burgt, I., Buske, A., Gillessen-Kaesbach, G., et al. (2008). Mutation and phenotypic spectrum in patients with cardio-facio-cutaneous and Costello syndrome. Clin. Genet. 73, 62–70. doi: 10.1111/j.1399-0004.2007.00931.x

PubMed Abstract | Crossref Full Text | Google Scholar

Sharland, M., Burch, M., McKenna, W. M., and Paton, M. A. (1992). A clinical study of Noonan syndrome. Arch. Dis. Child. 67, 178–183. doi: 10.1136/adc.67.2.178

PubMed Abstract | Crossref Full Text | Google Scholar

Tajan, M., Paccoud, R., Branka, S., Edouard, T., and Yart, A. (2018). The RASopathy family: consequences of germline activation of the RAS/MAPK pathway. Endocr. Rev. 39, 676–700. doi: 10.1210/er.2017-00232

PubMed Abstract | Crossref Full Text | Google Scholar

Tartaglia, M., Zampino, G., and Gelb, B. D. (2010). Noonan syndrome: clinical aspects and molecular pathogenesis. Mol. Syndromol. 1, 2–26. doi: 10.1159/000276766

Crossref Full Text | Google Scholar

Welch, B. L. (1947). The generalization of ‘stundent’s’ problem when several different population variances are involved. Biometrika 34, 28–35. doi: 10.1093/biomet/34.1-2.28

PubMed Abstract | Crossref Full Text | Google Scholar

Wendt, H. W. (1972). Dealing with a common problem in social science: a simplified rank-biserial coefficient of correlation based on the U statistic. Eur. J. Soc. Psychol. 2, 463–465. doi: 10.1002/ejsp.2420020412

Crossref Full Text | Google Scholar

Yoon, S., and Seger, R. (2006). The extracellular signal-regulated kinase: multiple substrates regulate diverse cellular functions. Growth Factors 24, 21–44. doi: 10.1080/02699050500284218

Crossref Full Text | Google Scholar

Zenker, M., and Kutsche, K. (2016). RASopathies. Med. Genet. 28, 15–38. doi: 10.1007/s11825-016-0080-8

Crossref Full Text | Google Scholar

Glossary