Molecular and Functional Characterization of Different BrainSphere Models for Use in Neurotoxicity Testing on Microelectrode Arrays

The currently accepted methods for neurotoxicity (NT) testing rely on animal studies. However, high costs and low testing throughput hinder their application for large numbers of chemicals. To overcome these limitations, in vitro methods are currently being developed based on human-induced pluripotent stem cells (hiPSC) that allow higher testing throughput at lower costs. We applied six different protocols to generate 3D BrainSphere models for acute NT evaluation. These include three different media for 2D neural induction and two media for subsequent 3D differentiation resulting in self-organized, organotypic neuron/astrocyte microtissues. All induction protocols yielded nearly 100% NESTIN-positive hiPSC-derived neural progenitor cells (hiNPCs), though with different gene expression profiles concerning regional patterning. Moreover, gene expression and immunocytochemistry analyses revealed that the choice of media determines neural differentiation patterns. On the functional level, BrainSpheres exhibited different levels of electrical activity on microelectrode arrays (MEA). Spike sorting allowed BrainSphere functional characterization with the mixed cultures consisting of GABAergic, glutamatergic, dopaminergic, serotonergic, and cholinergic neurons. A test method for acute NT testing, the human multi-neurotransmitter receptor (hMNR) assay, was proposed to apply such MEA-based spike sorting. These models are promising tools not only in toxicology but also for drug development and disease modeling.


Introduction
The currently accepted methods for neurotoxicity (NT) testing rely on animal studies. These are defined in the OECD test guidelines (TG) 418, 419, and 424 [1 -3]. The drawbacks of these TG studies are their resource intensities regarding money and time, their high variability, their lack of mechanistic understanding, and their uncertainties due to species differences [4,5]. Such species specificities between humans and rodents decrease confidence in TG neurotoxicity testing as the human brain holds unique features [6,7]. The human uniqueness of the brain is also reflected in the high attrition rates in central nervous system-related drug development when moving from preclinical research to clinical drug applications in humans [8,9]. Species-specific differences in pharmacokinetics and pharmacodynamics are one reason for this poor human prognosis [10]. Moreover, there is a lack of human-relevant neural disease models that also contributes to the high drug attrition rate [11].
One of the main arguments for a paradigm change in human health risk assessment, however, is the low testing throughput that has been leading to a large number of untested the dual SMAD inhibitors SB-431542 and LDN-193189 and supplemented with 20 ng/mL recombinant human basic FGF (#233-FB, R&D Systems, Wiesbaden, Germany) and 10 µM Y-27632 (only for the first 24 h after passaging). On day 21, hiNPCs were singularized with Accutase and frozen in neural progenitor medium containing 10% dimethyl sulfoxide (DMSO, #A994.1, Carl-Roth, Karlsruhe, Germany) and 10 µM Y-27632.

Thawing of hiNPCs
The vials of hiNPCs were quickly thawed in the palms of the hand and each vial containing 4 × 10 6 cells was directly diluted in 10 mL of the respective neural progenitor medium with 10 µM Y-27632 (#HB2297, Hello Bio, Bristol, UK). After centrifugation (300× g, 5 min), the cell pellet was resuspended in the respective neural progenitor medium with 10 µM Y-27632 (#HB2297, Hello Bio, Bristol, UK). The cells of one frozen cryovial were divided into three wells of a coated 6-well plate. The medium was replaced daily without the addition of Y-27632 (2D-NIM and GNEIB protocol) or with 10 µM Y-27632 (Stemdiff protocol).

Neural Differentiation on Microelectrode Arrays (MEA)
To access the neuronal electrical activity, the BrainSpheres were plated on 96-well multielectrode arrays (MEA, #M768-tMEA-96B, Axion Biosystems, Atlanta, GA, USA) after 3 weeks, 2 weeks, 1 week, or without 3D differentiation under constant shaking in CINDA+ or Electro differentiation medium. The MEA was coated with specific matrices for each differently generated and differentiated BrainSphere (see Supplementary Table S1). After coating with the respective matrix, one BrainSphere was placed in the middle of each well of the MEA, except for BrainSpheres generated with the GNEIB protocol and differentiated in Electro differentiation medium. Here, three BrainSpheres were placed per well as they were smaller in size. The BrainSpheres were fed twice per week by replacing half of the differentiation medium. Supplementary Figure S5B shows representative placements of the BrainSpheres on the MEAs after 3 weeks of 3D-differentiation and a further 3 weeks differentiated on the MEAs. The neuronal electrical activity was recorded twice per week, and BrainSpheres were acutely exposed to L-glutamate (50 µM), DL-2-Amino-5-phosphonovaleric acid (AP5, 50 µM), and NBQX disodium salt (NBQX, 50 µM), bicuculline (3 and 10 µM), picrotoxin (5, 10 and 20 µM), haloperidol (1 and 10 µM), carbaryl (5,10,20,50, and 100 µM), and buspirone hydrochloride (5,10,20,50, and 100 µM) after 2-6 weeks on the MEA. For the substance testing, the BrainSpheres were first exposed to the neurotransmitters glutamate (50 µM) or GABA (10 µM) before the antagonists AP5 (50 µM) and NBQX (50 µM) or bicuculline (10 µM) were added. The substances were removed with a complete exchange of the medium and the neuronal networks were allowed to recover for 2 to 3 h. After another baseline recording, the test compounds trimethyltin chloride (TMT) and emamectin were consecutively added until the final concentration was reached. Detailed information such as CAS registry numbers (CASRN), suppliers, and solvents is available in Supplementary Table S2. The data for subtype characterization were derived from 3 different MEA plates with 8 wells per condition and 8 electrodes per well, resulting in 192 electrodes per condition. The data of the substance testing experiment with TMT, emamectin, and quinpirole were derived from one MEA plate with 8 wells per condition and subtype specification, resulting in 64 electrodes.

Recording and Data Analysis of MEA Neuronal Electrical Activity
Extracellular recording of the neuronal electrical activity was performed twice a week for 15 min (baseline and each tested compound concentration) at 37 • C and 5% CO 2 after the cells were allowed to equilibrate for 15 min in the Axion Maestro Pro system (Axion Biosystems, Atlanta, GA, USA). Data recording was operated by the Axion Integrated Studios (AxIS) navigator software (version 3.1.2, Axion Biosystems, Atlanta, GA, USA) with a sampling frequency of 12.5 kHz and a digital band-pass filter of 200-3000 Hz. Subsequent spike detection was performed using the method "adaptive threshold crossing" with a threshold of 6 root mean square (rms) noise on each electrode and a pre-and postspike duration of 0.84 ms and 2.16 ms, respectively, and an electrode was defined as active with at least 2 spikes per min. Quantification of general electrical activity and neuronal network activity was performed with the Neural Metric Tool software (version 3.1.7, Axion Biosystems, Atlanta, GA, USA). If several spikes occur one after the other on the same electrode and meet certain criteria, they are referred to as a burst. For burst detection, the method "Inter-spike interval (ISI) threshold" was used with a minimum of five contributing spikes and a maximum ISI of 100 ms. If bursts occur simultaneously on several electrodes within a well, this is defined as a network burst. Network bursts were identified using the algorithm "envelope" with a threshold factor of 1.25, a minimal inter-burst interval (IBI) of 100 ms, at least 35% participating electrodes, and 75% burst inclusion. Parameters for neuronal activity (percentage of active electrodes and weighted mean firing rate (wMFR)) as well as for network maturation and synchronicity (burst frequency and network burst frequency) were analyzed.

Spike Sorting
For spike sorting, the AxIS generated .spk files of the baseline measurement and each corresponding treatment concentration were concatenated and converted into .nex files with a MATLAB (R2021b, R2022b, MathWorks, Natick, MA, USA) script. The generated .nex files were sorted with the neural spike sorting software Offline Sorter (OFS, version 4.4, Plexon, Dallas, TX, USA) using the automatic clustering T-Distribution EM method with 10 degrees of freedom (D.O.F) and 20 initial number of units. The sorted units were Cells 2023, 12, 1270 7 of 37 exported as per-unit and per-waveform data giving information about the number of spikes per unit per baseline and substance concentration. For analysis, only units with at least 2 spikes/min during the baseline measurement were analyzed and they were considered responding units when the fold change to the baseline was at least ±0.25 (increase or decrease).

Cytotoxicity Assessment
Cytotoxicity was assessed by measuring the lactate dehydrogenase (LDH) release from cells with damaged membranes using the CytoTox-ONE Homogeneous Membrane Integrity assay (#7891, Promega, Madison, WI, USA). For this, parallel to the substance testing on the MEA, one Brainphere was placed in the middle of each well of a coated 96-well plate (for coating see Supplementary Table S1) and differentiated in CINDA+ for 4 weeks. After acute treatment with the respective substance (15 min at 37 • C), 50 µL medium from each well was transferred to a new 96-well plate and 50 µL CytoTox-ONE reagent was added. As lysis control, neurospheres were treated with 10% Triton-X 100. The medium without spheres was used to correct for background fluorescence. Fluorescence was detected with the Tecan infinite M200 Pro reader (ex: 540 nm; em: 590 nm).

Quantitative Polymerase Chain Reaction (qPCR)
For gene expression analyses, samples were collected at different time points indicated in Figure 1. Messenger RNA (mRNA) was isolated using the Rneasy Mini Kit (#74104, Qiagen, Hilden, Germany) and transcribed into cDNA using the QuantiTect Reverse Transcription Kit (#205311, Qiagen, Hilden, Germany). qPCR was performed in the Rotor-Gene Q Cycler (#9001560, Qiagen, Hilden, Germany) using the QuantiFast SYBR Green PCR Kit (#204056, Qiagen, Hilden, Germany). For quantification, standard curves of all examined genes were generated for calculating copy numbers (CN) as described in Walter et al.'s work, 2019 [73]. All steps were performed according to the manufacturer's instructions and CN of the gene of interest were normalized to β-actin expression. Primer sequences are listed in Supplementary Table S3. The data were derived from three-four independent experiments, each performed with a different batch of hiNPCs.
Rotor-Gene Q Cycler (#9001560, Qiagen, Hilden, Germany) using the QuantiFast SYBR Green PCR Kit (#204056, Qiagen, Hilden, Germany). For quantification, standard curves of all examined genes were generated for calculating copy numbers (CN) as described in Walter et al.'s work, 2019 [73]. All steps were performed according to the manufacturer's instructions and CN of the gene of interest were normalized to β-actin expression. Primer sequences are listed in Supplementary Table S3. The data were derived from three-four independent experiments, each performed with a different batch of hiNPCs.  After thawing, the hiNPCs were cultivated in 2D to recover and proliferate. On day 4, they were transferred to a shaking incubator, formed BrainSpheres (3D), and were differentiated in CINDA+ or Electro medium for 1 week (3D-1w), 2 weeks (3D-2w), or 3 weeks (3D-3w). Characterization of gene expression (PCR), protein expression (ICC), and electrical activity (MEA) was performed at the indicated time points. ICC, immunocytochemistry; MEA, microelectrode array. Created with biorender.com.

Immunocytochemistry (ICC) of Adherent hiNPCs
On day 4 after thawing, hiNPCs were singularized with Accutase (10 min at 37 • C and 5% CO 2 ), and 20,000 cells were transferred into each chamber of a coated 8-chamber slide (#354118, BD Falcon, Franklin Lakes, NJ, USA) for another 3 days of cultivation before fixing with 4% paraformaldehyde (  After washing twice with DPBS, the BrainSpheres were transferred to a microscopy slide, Aqua-Poly/Mount (#18606-20, Polysciences Inc., Warrington, PA, USA) was added, and they were covered with a cover glass. Image acquisition was performed with the confocal laser scanning system LSM 710 (Zeiss, Oberkochen, Germany) at the Center for Advanced Imaging (CAi) of the Heinrich-Heine-University in Düsseldorf.

Statistical Analysis
All statistical analyses were performed with GraphPad Prism (version 8.4.3, Boston, MA, USA). The MEA data were examined for significant differences between the SDs using the Brown-Forsythe test. If the SDs were not significantly different, one-way ANOVA followed by the post hoc Dunnett test were applied. If the SDs were significantly different, Brown-Forsythe and Welch ANOVA followed by the post hoc Games-Howell test were performed. Flow cytometry data were analyzed using a two-way-ANOVA followed by Tukey test for correction of multiple comparisons. qPCR data were analyzed in two ways. For media comparisons, an unpaired t-test with Welch's correction was used. For comparison to the first time point (2D or 3D), one-way ANOVA and the post hoc test Dunnett were used. The calculated p-values for each comparison can be found in Supplementary Table S4.

All Three Neural Induction Protocols Successfully Induce hiPSCs into the Neural Lineage
The hiPSC line IMR-90 (clone 4, WiCell, Madison, WI, USA) was quality-controlled and banked as described in Tigges et al.'s work 2021 [69]. From a full quality-controlled master cell bank (MCB), a quality-controlled working cell bank (WCB) was prepared to ensure equal starting material for all experiments. The quality controls included the microscopic assessment of colony and cell morphology, mycoplasma detection, STR genotyping, caryotype analysis, protein expression of stem cell markers, viability assessment, gene expression analysis (only MCB), pluripotency assays (only MCB), and post-thaw recovery analysis. The data are shown in Tigges et al.'s work 2021 [69]. Three different neural induction protocols were compared for their capacity to generate human-induced neural progenitor cells (hiNPCs). Two of them were modified from previously published protocols (2D-NIM, [56]; GNEIB, [46,48]) and a third one was commercially available (Stemdiff). The hiPS cell line IMR90 was neurally induced as described in Figure 1A and characteristic molecular marker expression was analyzed using flow cytometry on day 12 (all three protocols) and day 21 (2D-NIM and GNEIB) of the protocol ( Figure 1B). The stemness marker octamer-binding transcription factor 3/4 (OCT3/4) [40] was downregulated in all protocols to under 4% positive cells. After 12 days of neural induction, the neural progenitor cell marker NESTIN [74] was expressed in 97% (on day 21: 99%), 94% (on day 21: 98%), and 99% of hiNPCs generated with the 2D-NIM, GNEIB, and Stemdiff protocols, respectively. Expression of the neural progenitor marker PAX6 [75] was unequally in the three protocols, i.e., in 93% and 98% of the cells generated with the GNEIB (day 12) and the Stemdiff protocol, respectively, but only in 58% of cells generated with the 2D-NIM (day 12) protocol. Additionally, the number of cells expressing PAX6 decreased on day 21 for 2D-NIM (15%) and GNEIB (67%). The proliferation marker KI-67 was expressed in 86% (2D-NIM), 87% (GNEIB), and 95% (Stemdiff) of the cells on day 12 and decreased to 66% (2D-NIM) and 67% (GNEIB) cells on day 21 of neural induction. After neural induction and FACS evaluation, the hiNPCs were frozen in liquid nitrogen to start each protocol with the same passage number and reduce inter-experimental variability. To ensure that the freezing process did not change molecular marker expression, we confirmed NESTIN and KI-67 expression after thawing via immunocytochemistry (Supplementary Figure S1). For forming BrainSpheres in 3D, the adherent hiNPCs were cultivated under shaking conditions in 6-well plates according to Honegger et al. 1979 andPamies et al., 2017 [57,76] using the two differentiation media, i.e., CINDA+ and Electro. These media were chosen for optimization of BrainSphere electrical activity. Their potential to support oligodendrocyte differentiation will be investigated in future studies. Spheres were characterized with regard to gene and protein expression, as well as neuronal network activity on MEAs in the proliferating state and after 1, 2 and 3 weeks of 3D differentiation ( Figure 1C).

BrainSpheres Differ in Neural Marker Gene Expression Depending on the Applied Protocol
For characterizing BrainSpheres' brain region specificity, developmental stage, neural subtypes, neurotransmitter processing and astrocyte maturation gene expression analyses were performed via qPCR at different time points as indicated in Figure 1. As markers for different brain regions, the genes FOXG1 (forkhead box G1, forebrain), LMX1A (LIM homeobox transcription factor 1 alpha, midbrain), EN1 (engrailed homeobox 1, midbrain), and HOXA2 (homeobox A2, hindbrain) were analyzed ( Figure 2, Supplementary Table S4) [43,77]. The forebrain marker FOXG1 was significantly higher expressed in BrainSpheres generated with the GNEIB protocol and almost not expressed in BrainSpheres generated with the 2D-NIM protocol. Although FOXG1 expression in Stemdiff-derived BrainSpheres was the highest, its induction was not statistically significant due to high standard deviations. The midbrain marker LMX1A was highest expressed in BrainSpheres generated with the GNEIB protocol; however, this was not statistically significant. The midbrain marker EN1 was higher expressed in BrainSpheres derived from the 2D-NIM and GNEIB protocols compared to BrainSpheres derived from the Stemdiff protocol, whose EN1 expression decreased significantly during differentiation. The hindbrain marker HOXA2 registered the highest expression in the early hiNPC stages of the 2D-NIM and GNEIB protocol, whereas the expression decreased significantly during three weeks of differentiation.
The genes SLC12A2 (solute carrier family 12 member 2) and SLC12A5 (solute carrier family 12 member 5) are encoding for the two ion transporters, i.e., Na + -K + -2Cl − cotransporter-1, (NKCC1) and K + -Cl − cotransporter (KCC2), respectively, that are key in the postnatal shift from excitatory to inhibitory GABAergic neurons caused by a reduced expression of SLC12A2 and an increased expression of SLC12A5 [78]. While SLC12A2 expression did not decrease significantly, SLC12A5 significantly increased in BrainSpheres derived from the 2D-NIM protocol; yet, with overall very low copy number expression. For comparison, mRNA expression from a fairly mature, 35 days differentiated neuron-glia mixed culture 2D network, i.e., the SynFire kit [72], was analyzed. SLC12A2 copy number expression was much lower and SLC12A5 copy number expression was much higher in the differentiated SynFire cells than in the BrainSpheres indicating prematurity of up to 3 weeks differentiated BrainSpheres independent of the applied protocol.
Expression of MAP2 (microtubule associated protein 2), which encodes for a dendritic protein [79], significantly increased during BrainSphere differentiation generated with the 2D-NIM and GNEIB protocols. The pre-synaptic marker SYN1 (synapsin 1) registered the highest expression in BrainSpheres generated with the 2D-NIM and the Stemdiff protocols and differentiated in CINDA+; however, this increase is not statistically significant as in 2D-NIM-BrainSpheres differentiated in Electro medium. The postsynaptic DLG4 (discs large MAGUK scaffold protein 4) expression was abundant in all BrainSpheres even without differentiation and reached similar or even higher expression levels (Stemdiff Electro) than in 35-day differentiated SynFire cells during differentiation.
For analyzing the potential for BrainSphere glial differentiation, the genes GFAP (glial fibrillary acidic protein), S100B (S100 calcium binding protein B) and AQP4 (aquaporin 4), which represent different stages of astrocyte maturation, were chosen [80]. GFAP was only expressed in BrainSpheres generated with the 2D-NIM protocol and differentiated in CINDA+ for at least 2 weeks. S100B was already expressed after one week of differentiation with significantly increased expression in BrainSpheres generated with the 2D-NIM (CINDA+) and the GNEIB protocol (both differentiation media), whereas both conditions differentiated in CINDA+ reached higher expression levels than the 35-day differentiated SynFire cells. AQP4, which denotes astrocytes with higher maturity, was hardly expressed in all BrainSpheres, yet highly expressed in SynFire cells.
Expression of genes encoding for the glutamate receptors GRIA1 (glutamate ionotropic receptor AMPA type subunit 1) and GRIN1 (glutamate ionotropic receptor NMDA type subunit 1) were only significantly increased during differentiation in BrainSpheres generated with the 2D-NIM protocol in both differentiation media; however, the GRIA1 expression in Stemdiff BrainSpheres exceeded the expression of 35-day differentiated SynFire cells. The expression of GABRA1 (gamma-aminobutyric acid type A receptor subunit alpha 1) increased during differentiation in BrainSpheres generated with the 2D-NIM and the Stemdiff protocols. The gene DRD2 encoding for the dopamine receptor D2 was not significantly expressed; however, it registered higher expression levels in BrainSpheres generated with the GNEIB protocol. The genes HTR1A (5-hydroxytryptamine receptor 1A) and CHRNA4 (cholinergic receptor nicotinic alpha 4 subunit), encoding for serotonin and choline receptors, respectively, were scarcely expressed in all conditions. All receptors were expressed in the SynFire cells.  In addition to neurotransmitter receptors, synthesis and transport of neurotransmitters are also crucial for neural functioning. The enzyme glutaminase (GLS) is necessary for glutamate synthesis [81] and the encoding gene was already expressed in the undifferentiated BrainSpheres and increased significantly only in GNEIB-BrainSpheres during differentiation in CINDA+, whereas the gene for vesicular glutamate transporter 1 (SLC17A7) registered the highest expression levels in Stemdiff BrainSpheres after 2 and 3 weeks of differentiation in CINDA+; however, this was not statistically significant. The expression of vesicular glutamate transporter 2 (SLC17A6) did not strongly increase during differentiation. Glutamate decarboxylase 1 (GAD1), which catalyzes a critical step in GABA synthesis [82], registered the highest expression levels in BrainSpheres generated with the 2D-NIM protocol and the expression significantly increased during maturation. The highest expression levels of tyrosine hydroxylase (TH) and tryptophan hydroxylase 1 (TPH1) were found in GNEIB-BrainSpheres differentiated in CINDA+. The genes CHAT (choline acetyltransferase) and ACHE (acetylcholinesterase) were hardly expressed in all Brain-Spheres but highly so in the SynFire cells. Additionally, SLC17A7, GAD1, TH, and TPH1 registered higher expression levels in some BrainSpheres than in the 35-day differentiated SynFire cells.

Neural Induction and Differentiation Protocols Determine the Potential of BrainSpheres to Differentiate into Astrocytes and Dopaminergic Neurons
Immunofluorescence stainings of BrainSpheres differentiated for 3 weeks revealed their differentiation potentials into astrocytes (S100B) and dopaminergic neurons (TH). BrainSpheres neurally induced with the 2D-NIM protocol differentiated into an abundance of S100B-positive cells that formed a dense layer underneath the TUBB3-positive neurons, regardless of the differentiation medium used. In contrast, BrainSpheres derived from the other two neural induction protocols generated only a few or no cells of the astrocytic lineage ( Figure 3A). Furthermore, BrainSpheres generated with the GNEIB protocol and differentiated in CINDA+ generated the most TH-positive dopaminergic neurons followed by BrainSpheres generated with the 2D-NIM protocol (CINDA+ and Electro) ( Figure 3B). The other BrainSphere conditions showed no or only a few TH-positive neurons. Furthermore, the BrainSpheres had different sphere and neuron morphologies. Especially of note, GNEIB-BrainSpheres differentiated in Electro medium showed a more inhomogeneous distribution of neurons.

Neural Induction and Differentiation Media Determine Neuronal Activity and Neural Network Function of BrainSpheres on MEAs
MEAs are powerful tools to evaluate the electrical activity of neuronal networks providing information on, e.g., the number of active electrodes, the rate of action potentials per electrode (weighted mean firing rate; wMFR), the clustering of spikes to bursts as a sign of neuronal maturation (burst frequency) and the neuronal network activity mirroring synchronous communication between different neurons within the neuronal network (network burst frequency). To evaluate and compare the functionality of the neural networks (NN), BrainSpheres were plated on MEAs either without or after 3D differentiation on a gyrical shaker for 1, 2, or 3 weeks, and the electrical activity was measured twice a week for 7 additional weeks. In general, the most active electrodes and the highest wMFR were observed after 3 weeks of 3D differentiation before plating BrainSpheres on the MEAs (Figure 4, Supplementary Figures S2-S4). Therefore, all of the following experiments were performed with this condition. The number of active electrodes is a valuable parameter as we observed that an inactive electrode does not necessarily indicate that neurons are not growing across the electrode; however, a non-electrically active cell may cover it (Supplementary Figure S5). The NN generated with the 2D-NIM induction protocol and 3D differentiated in CINDA+ medium showed the most active wells (83%) and electrodes (58%) with the highest wMFR of 6.86 Hz after 7 weeks on MEAs ( Figure 4, Table 1). BrainSpheres derived from the GNEIB protocol differentiated in CINDA+ had the second most active wells (83%), electrodes (31%), and wMFR (6.30 Hz), but also had the highest variance with a standard error of mean (SEM) of 1.25 Hz. Both NN showed a similar mean burst frequency (2D-NIM: 0.31 Hz, GNEIB: 0.36 Hz) and an increasing activity over the course of 7 weeks on the MEAs. In contrast, NN derived from BrainSpheres either neurally induced with the 2D-NIM or GNEIB protocol and differentiated in Electro medium, or neurally induced with the Stemdiff protocol and differentiated in CINDA+ medium showed a decrease in the wMFR and burst frequency over the 7 weeks starting from 4 weeks on the MEAs. The NN derived from BrainSpheres neurally induced with the GNEIB protocol and differentiated in CINDA+ had the highest network burst frequency with 0.15 Hz. The percentages of spikes contributing to a network burst (network burst percentage) were similar for all conditions that generated network bursts after 7 weeks on MEA, but highest for the NN generated with the Stemdiff neural induction protocol and subsequent differentiation in Electro medium (48.02%, Table 1). In general, neural induction of hiPSCs using the 2D-NIM and GNEIB protocols followed by BrainSphere differentiation in CINDA+ for 3 weeks generated the most functional NNs with the most active electrodes and highest wMFR, burst frequencies and network burst frequencies, indicating that these protocols are favorable for the generation of mature NN from BrainSpheres.    The culture conditions with the strongest responses to each treatment are highlighted in bold. wMFR, weighted mean firing rate; glu, glutamate; bic, bicuculline; ptx, picrotoxin; halo, haloperidol; bsp, buspirone hydrochloride; crb, carbaryl.

Neural Induction and Differentiation Media Determine BrainSpheres' Neuronal Subtype Differentiation
For characterization of the NN not only the general electrical activity is important but also the responses to specific pharmacological modulators which highly depends on the occurrence of different neuronal subtypes. Therefore, we addressed the presence of glutamatergic, GABAergic, dopaminergic, serotonergic, and cholinergic neurons in 3-week 3D differentiated BrainSpheres after 2 to 6 weeks on MEAs. Quantification of neuronal units' wMFR responding to pharmacological modulation was possible due to spike sorting with the software Offline Sorter ( Figure 5A). This allowed the identification of specific responses of individual neurons within the integrated neuronal activities of single electrodes. Neural units reacting to the modulation with a change of at least ±25% in comparison to the baseline measurement were defined as responding (glutamatergic, GABAergic, dopaminergic, serotonergic, or cholinergic) units. The percentage of responding units, their fold-change to the untreated baseline measurement (colored), and the comparison to the fold-changes of the unsorted (grey) units were analyzed.  towards GABA receptor antagonism ( Figure 6). The raw values of the wMFR after exposure to bicuculline an picrotoxin are shown in Supplementary Figures S7 and S8.  First, the glutamatergic response was measured by applying the neurotransmitter glutamate followed by the application of the two glutamate receptor antagonists AP5 (antagonizes NMDA receptors) and NBQX (blocks AMPA receptors) [83]. It is expected that glutamate increases the electrical activity as an excitatory neurotransmitter. Without spike sorting, glutamate-dependent increased electrical activity was only measured in NN derived from BrainSpheres differentiated with GNEIB/Electro media, whereas BrainSpheres derived from the other five protocols showed no clear positive response or even a decreased wMFR upon acute glutamate exposure ( Figure 5B). After spike sorting, glutamate-responsive units from BrainSpheres produced with the 2D-NIM and GNEIB protocols produced the highest response to glutamate. Half of glutamate-responding neuronal units in 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres also responded with a decline in the wMFR to subsequent glutamate receptor inhibition by AP5/NBQX. In contrast, wMFRs of all neuronal glutamate-responsive units from BrainSpheres produced with Stemdiff/CINDA+ and Stemdiff/Electro were inhibited by AP5/NBQX. However, the latter two protocols per se produced very few responding units with minor changes in the wMFR amplitude ( Figure 5B, Table 1). In addition to the fold changes shown in Figure 5B, the raw values of the wMFR are presented in the supplementary material (Supplementary Figure S6).
To characterize the GABAergic response, the GABA receptor antagonists bicuculline and picrotoxin were applied. Bicuculline binds to GABA A receptors, whereas picrotoxin targets GABA A and GABA C receptors [84,85]. It is expected that upon treatment with the GABA receptors antagonists, mature GABAergic neurons post the GABA switch respond with an increased activity, e.g., wMFR, while immature GABAergic neurons before the GABA switch respond with a decreased firing. Without spike sorting, none of the Brain-Sphere cultures derived from any of the six protocols produced changes in NN activity upon GABA receptor antagonism ( Figure 6). After spike sorting, neuronal units of all BrainSphere protocols were identified that contained increased or decreased wMFR upon bicuculine and picrotoxin treatment; yet, most of the responses showed considerable variation in magnitude of responses and thus lacked statistical significance. Only NN from BrainSpheres generated with 2D-NIM/CINDA+, 2D-NIM/Electro, Stemdiff/CINDA+ and Stemdiff/Electro significantly increased the wMFR upon bicuculine treatment, with NN from 2D-NIM/Electro generated BrainSpheres being the least sensitive with the earliest response at 10 µM. In response to picrotoxin only the BrainSpheres generated with 2D-NIM/CINDA+ (5 µM) significantly induced the wMFR of NN. Neuronal units responding with an increased wMFR towards the GABA receptor antagonists ranged between 16% (Stemdiff/CINDA+) and 62% (2D-NIM/Electro) of all active neurons. However, the 2D-NIM/CINDA+ media generated the highest absolute number of positively responding neuronal units ( Figure 6). In addition to the units reacting to the GABA antagonists with increased wMFR, some units responded with a decreased activity. This is in line with the low gene expression of SLC12A5, which encodes for the ion receptor KCC2 (Figure 2). KCC2 increases its expression after the postnatal GABA switch from excitatory to inhibitory GABAergic neurons. Therefore, without spike sorting, the pre-mature and the more mature GABA receptor containing units compensated each other and resulted in no visible response to the antagonists. BrainSpheres generated with 2D-NIM/CINDA+ media produced the highest absolute numbers of neuronal units responding with an inhibitory and decreased wMFR towards the GABA receptor antagonists. Regarding the Stemdiff neural induction differentiated in CINDA+, in relation to all firing units, 31% (picrotoxin) to 59% (bicuculline) of BrainSphere-derived units responded with a decrease in activity towards GABA receptor antagonism ( Figure 6). The raw values of the wMFR after exposure to bicuculline an picrotoxin are shown in Supplementary Figures S7 and S8.
Inhibitory dopaminergic D2 receptors were addressed by applying the antagonist haloperidol to the cultures, which should increase the wMFR [86,87]. Without spike sorting, haloperidol enhanced the wMFR in GNEIB/CINDA+ (Figure 7). After spike sorting, all protocols derived units responding with an increased and decreased wMFR after exposure to haloperidol, except for BrainSpheres from 2D-NIM/Electro media, which only resulted in units with decreased activity. Haloperidol decreased the wMFR of neuronal units in NN generated with 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres most effectively with significant wMFR reductions at 1 µM haloperidol in 57% and 65% of all neuronal units, respectively. In addition, NN generated with the GNEIB/CINDA+ protocol showed the strongest increased activity at 1 µM. However, most BrainSphere neurons differentiated with 2D-NIM/CINDA+ media seem to express inhibitory D2 receptors as 20% of all recorded neuronal units responded to 1 µM haloperidol with an increased activity (Figure 7). Rising haloperidol concentrations decrease the number of units responding with an increase in wMFR, while the number of units reacting with decreased wMFR enhances under the treatment. The raw wMFR values of the unsorted and the responding units after exposure to haloperidol are shown in Supplementary Figure S9 and Table 1.
Buspirone, which agonizes the inhibitory serotonin 5-HT1A receptor, was applied to investigate the presence of serotonergic responses in the NN [88]. Without spike sorting, buspirone did not cause significant changes in neuronal activity in any of the BrainSpheres. However, after sorting, it decreased the wMFR of NN generated with 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres most effectively with significant wMFR reductions at 5 µM buspirone. In BrainSpheres produced using these protocols, 26% and 52% of all neuronal units responded to the compound with a decreased activity, respectively. BrainSpheres differentiated in Electro medium overall produced fewer serotonergic neurons than CINDA+ differentiated spheres with all neural induction protocols. However, with 88% of all units, neurons derived from BrainSphere generated with 2D-NIM/Electro media exhibited the highest percentage of neurons responding to buspirone with a decreased activity (Figure 8, Table 1). All protocols generated lower percentages of units responding with an increased wMFR to buspirone compared to a decreased activity except for BrainSpheres generated with 2D-NIM/CINDA+, where 31% of units responded with an increase. In addition to the fold change, the raw values of the wMFR are shown in Supplementary Figure  Inhibitory dopaminergic D2 receptors were addressed by applying the antagonist haloperidol to the cultures, which should increase the wMFR [86,87]. Without spike sorting, haloperidol enhanced the wMFR in GNEIB/CINDA+ (Figure 7). After spike sorting, all protocols derived units responding with an increased and decreased wMFR after exposure to haloperidol, except for BrainSpheres from 2D-NIM/Electro media, which only resulted in units with decreased activity. Haloperidol decreased the wMFR of neuronal units in NN generated with 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres most effectively with significant wMFR reductions at 1 µM haloperidol in 57% and 65% of all neuronal units, respectively. In addition, NN generated with the GNEIB/CINDA+ protocol showed the strongest increased activity at 1 µM. However, most BrainSphere neurons Figure 6. Neuronal network characterization via acute pharmacological modulation to assess GABAergic response. BrainSpheres were 3-week 3D differentiated before plated on MEA and consecutively exposed to (A) bicuculline (bic) or (B) picrotoxin (ptx). Shown are the fold changes to the untreated baseline measurements of all units (unsorted) and the sorted responding units (colored). Data are represented as mean ± SEM of three independent MEA experiments with eight wells per condition (*: significant to unsorted, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). The numbers above the bars represent the number of units that responded accordingly. differentiated with 2D-NIM/CINDA+ media seem to express inhibitory D2 receptors as 20% of all recorded neuronal units responded to 1 µM haloperidol with an increased activity ( Figure 7). Rising haloperidol concentrations decrease the number of units responding with an increase in wMFR, while the number of units reacting with decreased wMFR enhances under the treatment. The raw wMFR values of the unsorted and the responding units after exposure to haloperidol are shown in Supplementary Figure S9 and Table 1. Buspirone, which agonizes the inhibitory serotonin 5-HT1A receptor, was applied to investigate the presence of serotonergic responses in the NN [88]. Without spike sorting, buspirone did not cause significant changes in neuronal activity in any of the Brain-Spheres. However, after sorting, it decreased the wMFR of NN generated with 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres most effectively with significant wMFR reductions at 5 µM buspirone. In BrainSpheres produced using these protocols, 26% and 52% of all neuronal units responded to the compound with a decreased activity, respectively. BrainSpheres differentiated in Electro medium overall produced fewer serotonergic neurons than CINDA+ differentiated spheres with all neural induction protocols. However, with 88% of all units, neurons derived from BrainSphere generated with 2D-NIM/Electro media exhibited the highest percentage of neurons responding to buspirone with a decreased activity ( Figure 8, Table 1). All protocols generated lower percentages of units responding with an increased wMFR to buspirone compared to a decreased activity except for BrainSpheres generated with 2D-NIM/CINDA+, where 31% of units responded with an increase. In addition to the fold change, the raw values of the wMFR are shown in Supplementary Figure S10.
Cholinergic signal transduction was modulated with the insecticide carbaryl, which inhibits the enzyme acetylcholinesterase (AChE) [89] and binds to nicotinic acetylcholine receptors (nAChR) [90]. While AChE blockage increases, binding to nAChR decreases cholinergic neuronal activity [91]. Without spike sorting, no significant change in wMFR was observed after exposure to 5 µM carbaryl (Figure 9). After spike sorting, carbaryl Figure 7. Neuronal network characterization via acute pharmacological modulation to assess dopaminergic responses. BrainSpheres were 3-week 3D differentiated before plated on MEAs and exposed to haloperidol (halo). Shown are the fold changes to the untreated baseline measurement of all units (unsorted, grey) and the responding units after sorting (colored). Data are represented as mean ± SEM of three independent MEA experiments with eight wells per condition (*: significant to unsorted, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). The numbers above the bars represent the total number of units that respond accordingly.
Cholinergic signal transduction was modulated with the insecticide carbaryl, which inhibits the enzyme acetylcholinesterase (AChE) [89] and binds to nicotinic acetylcholine receptors (nAChR) [90]. While AChE blockage increases, binding to nAChR decreases cholinergic neuronal activity [91]. Without spike sorting, no significant change in wMFR was observed after exposure to 5 µM carbaryl (Figure 9). After spike sorting, carbaryl decreased the wMFR in NN generated with 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres most effectively with significant wMFR reductions at 5 µM carbaryl and 22-30% of all neuronal units responding to the compound. Both Stemdiff/CINDA+ and Stemdiff/Electro protocols produced BrainSpheres with the lowest number of (11 and 13%, respectively) and least sensitive (50 and 10 µM carbaryl, respectively) responding neuronal units ( Figure 9A, Table 1). NN derived from BrainSpheres generated with the 2D-NIM/CINDA+ and GNEIB/CINDA+ protocols also exhibited the highest absolute number of neuronal units reacting with the highest increased activity in wMFR to carbaryl ( Figure 9B). In addition, 2D-NIM BrainSpheres were overall the most sensitive to wMFR modulation in both directions, significantly responding at 5 µM (Figure 9). Interestingly, rising carbaryl concentrations increased the number of neuronal units that responded with a decrease in wMFR, while the number of units reacting with enhanced wMFR decreased under the treatment. In addition to the fold change, the raw values of the wMFR are shown in Supplementary Figure S11. 3.6. Set-Up of a New NAM for Acute Neurotoxicity Testing Using MEAs and Spike Sorting, the Human Multi-Neurotransmitter Receptor (hMNR) Assay With the well-characterized BrainSpheres, we propose the set-up of a test method as a NAM for acute neurotoxicity testing using MEAs and spike sorting. While general MEA activity can provide an overview of the general changes in NN activity, its resolution is not high enough to understand individual neuronal responses. Spike sorting in combination with neuronal subtype-specific model compounds seems to be a valuable solution for neuronal subtype identification in BrainSpheres on MEAs. To study if this system is suitable for acute neurotoxicity assessment, we set up a standard operating procedure combining neuronal unit identification with consecutive compound testing ( Figure 10A). With this setup, as a proof-of-concept, we measured the effects of two compounds, i.e., TMT, which enhances gluatamate release, and emamectin, a GABA-receptor agonist [92,93], on glutamatergic and GABAergic neuronal units in differentiated BrainSpheres ( Figure 10B). These compounds are the first substances of a chemical training set for the test method and were selected from the mode-of-action analyses in Masjosthusmann et al.'s work 2018 [22]. BrainSpheres generated with the 2D-NIM/CINDA+ protocol were used due to the resulting higher number of active electrodes and lower variance in comparison to the other protocols. Neural units reacting to neurotransmitter and antagonist with a change of at least ±25% in comparison to the baseline measurement were defined as responding glutamatergic or GABAergic units (Supplementary Figure S12A). decreased the wMFR in NN generated with 2D-NIM/CINDA+ and GNEIB/CINDA+ BrainSpheres most effectively with significant wMFR reductions at 5 µM carbaryl and 22-30% of all neuronal units responding to the compound. Both Stemdiff/CINDA+ and Stemdiff/Electro protocols produced BrainSpheres with the lowest number of (11 and 13%, respectively) and least sensitive (50 and 10 µM carbaryl, respectively) responding neuronal units ( Figure 9A, Table 1). NN derived from BrainSpheres generated with the 2D-NIM/CINDA+ and GNEIB/CINDA+ protocols also exhibited the highest absolute number of neuronal units reacting with the highest increased activity in wMFR to carbaryl ( Figure 9B). In addition, 2D-NIM BrainSpheres were overall the most sensitive to wMFR modulation in both directions, significantly responding at 5 µM ( Figure 9). Interestingly, rising carbaryl concentrations increased the number of neuronal units that responded with a decrease in wMFR, while the number of units reacting with enhanced wMFR decreased under the treatment. In addition to the fold change, the raw values of the wMFR are shown in Supplementary Figure S11. EER REVIEW 25 of 38 ** p ≤ 0.01, *** p ≤ 0.001). The numbers above the bars represent the total number of units that respond accordingly. Figure 9. Neuronal network characterization via acute pharmacological modulation to assess cholinergic responses. BrainSpheres were 3-week 3D differentiated before plated on MEA and exposed to carbaryl (crb). (A) Decreased and (B) increased responses after exposure to crb were detected. Shown are the fold changes to the untreated baseline measurement of all units (unsorted, grey) and the responding units after sorting (colored). Data are represented as mean ± SEM of three independent MEA experiments with eight wells per condition (*: significant to unsorted, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). The numbers above the bars represent the number of units that respond accordingly.

Set-Up of a New NAM for Acute Neurotoxicity Testing Using MEAs and Spike Sorting, the Human Multi-Neurotransmitter Receptor (hMNR) Assay
With the well-characterized BrainSpheres, we propose the set-up of a test method as a NAM for acute neurotoxicity testing using MEAs and spike sorting. While general MEA activity can provide an overview of the general changes in NN activity, its resolution is not high enough to understand individual neuronal responses. Spike sorting in combination with neuronal subtype-specific model compounds seems to be a valuable solution for neuronal subtype identification in BrainSpheres on MEAs. To study if this system is suitable for acute neurotoxicity assessment, we set up a standard operating procedure combining neuronal unit identification with consecutive compound testing ( Figure 10A). With this setup, as a proof-of-concept, we measured the effects of two compounds, i.e., TMT, Figure 9. Neuronal network characterization via acute pharmacological modulation to assess cholinergic responses. BrainSpheres were 3-week 3D differentiated before plated on MEA and exposed to carbaryl (crb). (A) Decreased and (B) increased responses after exposure to crb were detected. Shown are the fold changes to the untreated baseline measurement of all units (unsorted, grey) and the responding units after sorting (colored). Data are represented as mean ± SEM of three independent MEA experiments with eight wells per condition (*: significant to unsorted, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). The numbers above the bars represent the number of units that respond accordingly.
After exposure to TMT and respective spike sorting, TMT caused an increased wMFR in glutamatergic units in the sub-micromolar range with a decreasing effect starting at 2.22 µM TMT ( Figure 10B). This expected increasing effect in wMFR due to enhanced glutamate release by TMT [93] was not observed in the unsorted and GABAergic units. Additionally, the highest TMT concentration (20 µM) also decreased the wMFR of the unsorted and the GABAergic units. Treatment with the GABA A and GABA C receptor agonist emamectin [92] decreased the wMFR of all unsorted and sorted units in a concentrationdependent manner. However, the effect was strongest in sorted GABAergic units with a reduction to a fold change of 0.34 at 30 nM ( Figure 10B). That the emamectin effects can even be observed in the unsorted units can be explained by the high abundance of GABAergic units in these particular BrainSpheres, with 80 to 95% of the units reacting to GABA and bicuculline.
Unspecific cytotoxic effects of the two test compounds were excluded by measuring LDH release (Supplementary Figure S12B).

Discussion
In recent years, industry, regulators, and academia have agreed on the need for NAMs to test chemicals with higher throughput, lower costs, and better predictivity for humans [13,15]. For this task, human cell systems designed for a specific purpose should Figure 10. Test method set-up for the NAM assessing neuronal subtype-specific acute neurotoxicity using MEA and spike sorting. (A) Schematic workflow for acute neurotoxicity testing. After the first baseline recording (baseline 1), first, the indicated neurotransmitter (glutamate or GABA), followed by the corresponding antagonist (AP5/NBQX or bicuculline), was applied. The pharmacological modulators were removed via a complete washout and the neuronal networks were allowed to recover for 2 to 3 h. After a second baseline measurement (baseline 2), the test substances TMT and emamectin were applied in increasing concentrations. (B) 2D-NIM BrainSpheres were 3 weeks 3D-differentiated in CINDA+. After 4 subsequent weeks of differentiation on the MEA, they were used for the described acute neurotoxicity testing. Shown are the fold changes of the wMFR to the respective untreated baseline measurements of all units (unsorted) or pre-sorted responding units (colored). The horizontal bar under each graph indicates the percentages of units responding to the modulation with neurotransmitter and antagonist and were thus defined as glutamatergic or GABAergic units. Data are represented as mean ± SEM (*: significant compared to the baseline, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). Created with biorender.com.

Discussion
In recent years, industry, regulators, and academia have agreed on the need for NAMs to test chemicals with higher throughput, lower costs, and better predictivity for humans [13,15]. For this task, human cell systems designed for a specific purpose should preferably be used and combined with other in vitro and in silico methods to cover multiple endpoints [94,95]. Human in vitro models for acute neurotoxicity testing that examine neurotransmission mainly refer to effects of total neural networks by measuring spike-, burst-, and network-related parameters. While these parameters provide valuable information, they do not necessarily account for a large variety of neurotoxic MoA [31,37,[96][97][98][99][100]. The low granularity of the classical MEA evaluation by studying integrated signals over single electrodes is accompanied by a high dependence on the NN composition. As hiPSCderived NN are especially fairly variable [43], high uncertainty might be involved in using spontaneously formed NN from hiPSC for in vitro neurotoxicity studies. Therefore, we characterized mixed culture BrainSpheres [57] for setting up a multiplexed test method for acute neurotoxicity evaluation with the goal of adding multiple neurotoxicity MoA to the established parameters measured with MEAs. First, we characterized six different BrainSphere models resulting from three adherent neural induction protocols combined with two different media for subsequent differentiation. All neural induction protocols showed low variability and high efficiency by resulting in at least 97% cells expressing the neural progenitor marker NESTIN [45,74]. However, hiNPC generated with the 2D-NIM, yet not with the GNEIB protocol, contained less KI-67 positive cells on day 21 of neural induction. This might be the consequence of asymmetrical cell division into proliferative and non-proliferative daughter cells [101]. PAX6 controls various processes regarding the neuroectodermal fate in a concentration-dependent manner and, if absent, leads to asymmetric cell division and thus to neurogenesis [75,102,103]. Hence, the low PAX6 expression in 2D-NIM hiNPCs could explain the decrease in proliferating cells and indicate a more developed state in comparison to the other two protocols.
After the successful neural induction, hiNPCs were frozen in liquid nitrogen. This allows each subsequent experiment to be performed with the same hiNPC passage number, reduces variability and saves time and money. We confirmed that this additional step does not alter the expression of NESTIN and KI-67.
Gene expression data at various differentiation times showed that BrainSpheres produced in different media differ in genes referring to distinct brain regions, synapse formation, astrocyte differentiation, receptors, and neuronal subtypes. Interestingly, the neural induction media influence gene expression more strongly than the differentiation media. Moreover, gene expression comparison to the SynFire neural cells, which are cell ratiocontrolled, pre-differentiated excitatory and inhibitory neurons and astrocytes forming functional and highly synchronous neural networks over 35 days in culture, revealed that the BrainSpheres are still fairly immature after 3 weeks of differentiation (e.g., comparing MAP2, SYN1, DLG4, AQP4 expression). This is supported by the low expression of the K + -Cl − -co-transporter (KCC2, SLC12A5) we observed in BrainSpheres compared to the SynFire cells. The increase in KCC2 expression together with the decrease in expression of the Na + -K + -2Cl − -co-transporter (NKCC1, SLC21A2) initiates the postnatal switch from excitatory to inhibitory GABA signaling [78]. The presence of mature and immature GABAergic neurons was also supported by the MEA measurements after exposure to the GABA antagonists bicuculline and picrotoxin, which resulted in both increases and decreases in wMFR. We observed very low expression of genes encoding for serotonin receptor (HTR1A), choline receptor (CHRNA4), and choline synthesis enzyme (CHAT); however, the MEA analysis revealed functional receptors. Previous studies showed that some neuronal subtype-specific markers were only expressed in mature neurons, which can take up to 16 weeks to achieve with hiPSC-derived NN [49,104]. Therefore, such gene expression data in mixed cultures have to be regarded with caution as gene expression measured is an integration of cellular expression and cell abundance in the cultures. Hence, protein analyses using immunocytochemistry combined with functional studies, e.g., using MEAs, are needed for proper test system characterization.
Immunocytochemical stainings for different neuronal and astrocytic markers revealed different potentials of the distinct neural induction/differentiation protocols for differentiation into S100B-positive cells of the astrocyte lineage, with the 2D-NIM/CINDA+ and GNEIB/CINDA+ protocols being the most effective, while GNEIB/Electro and Stemdiff/Electro BrainSpheres only differentiated into very few or no S100B positive cells. These data were supported by the gene expression analyses. Astrocytes are essential for NN maturation and function since they play an important role in synaptogenesis, neuronal survival and outgrowth, phagocytosis, and NT uptake from the synaptic cleft [105,106]. However, spatiotemporal astrocyte marker expression has to be considered when analyzing astrocytes in vitro. Data from human in vivo investigations reveal that GFAP, a marker most commonly used in in vitro studies, should not be used as a general and sole astrocytic marker [80]. Therefore, immunocytochemical analyses supported by qPCR results using a panel of astrocytic lineage markers (e.g., GFAP, S100B, and AQP4) seem reasonable. Astrocyte presence is also important concerning the effects of chemicals, as neurons and astrocytes might react differently to toxic substances [107][108][109]. Moreover, astrocytes might be the mediators of neuronal toxicity [110,111] or even neuroprotective [112,113]. Hence, astrocytes' presence in mixed co-cultures is thought to enhance the applicability domain compared to pure cultures.
Not only do astrocytes and neurons respond differently to certain chemicals, but toxic effects on individual neuronal subtypes are also frequent causes of neurotoxicity [22,59,114]. Prominent examples are the pesticide rotenone, which acts specifically on dopaminergic neurons [115], or the acetylcholinesterase inhibitor parathion targeting cholinergic neurons [116]. In that regard, it might be of interest to generate fit-for-purpose BrainSpheres enabling the study of distinct neuronal subtypes depending on the scientific or regulatory question. Here, for example, ICC stainings revealed that BrainSpheres generated with the 2D-NIM (both differentiation media) or the GNEIB/CINDA+ protocols generated the highest numbers of TH-positive dopaminergic neurons, while in Stemdiff BrainSpheres, they were much fewer. However, besides cell type-specific marker expression, for a physiologically relevant neural test systems, the formation of a functional neuronal network and adequate responses to model compounds have to be demonstrated.
In addition to the mRNA and protein expression data, BrainSpheres were characterized for their performance on MEAs. Similar to the expression analyses, the applied induction and differentiation protocols determined the BrainSphere's activities on the MEAs. BrainSpheres induced in 2D-NIM and GNEIB media differentiated in CINDA+ showed the most active electrodes, the highest wMFRs, and the highest burst frequencies.
The wMFR and bursting behavior depend on various factors, predominantly the presence of astrocytes and the ratio of excitatory to inhibitory neurons [31,32,34]. According to the expression data, these two protocols express the highest levels of the astrocytic lineage marker S100B. Hence, astrocyte presence might be responsible for the abundant firing activity of these BrainSpheres. In contrast to BrainSpheres neurally induced with the 2D-NIM and GNEIB media, BrainSphers generated with Stemdiff medium displayed higher electrical activity when subsequently differentiated with the Electro medium. This indicates that the combination of neural induction and differentiation media highly influences NN functionality. To date, only the influence of different neural induction protocols or various differentiation conditions have been analyzed; however, a combination of both has not been examined so far [53,56,117,118].
While spike, burst, and network parameters provide important information on general network function, their level of granularity is not particularly high. Therefore, we applied the method of spike sorting to the MEA data [119]. Spike sorting enables the identification of single active neuronal 'units' within the signal of one MEA electrode using curve progression analyses. These units can be evaluated individually and hence quantified across multiple electrodes. We challenged the NN with model compounds targeting glutamate, GABA, dopamine, serotonin and nACh receptors, as well as acetylcholine esterase to identify different neuronal subtype signalling.
Without spike sorting, glutamate did not significantly increase wMFR signals on MEAs. After spike sorting, neuronal units were identified that increased or decreased their activity. Neuronal response of increased activity was expected for an excitatory neurotransmitter. The opposite effect might be attributed to presynaptic metabotropic glutamate receptors (mGluR) that can act as a negative feedback loop and inhibit glutamate release [120][121][122][123][124] possibly in early developing neurons as a counterpart to immature excitatory GABAergic neurons [125]. Such opposite glutamate effects were observed earlier in other neural in vitro models [32,56].
Spike sorting of MEA activity after exposure to the GABA receptor antagonists bicuculline and picrotoxin revealed that all established NN contain both excitatory and inhibitory GABAergic neurons. This was not visible in the whole electrode recordings, since the two opposite reactions cancel each other out. Excitatory action of GABA is a physiological response before the GABA switch [78]. Hence, all NN also contain immature neurons that precede the GABA switch in addition to inhibitory GABAergic neurons. GABA A and GABA C receptors seem to mature differently as picrotoxin, which binds to GABA A and GABA C receptors [85], increases the increase/decrease ratios of neuronal units compared to bicuculline, which only interacts with GABA A receptors [84], suggesting higher maturation states of GABA A receptors. The gene expression for the two ion transporters NKCC1 and KCC2, which marks the switch from pre-mature excitatory to mature inhibitory GABAergic neurons, supports these functional data.
All BrainSphere models also responded to haloperidol, buspirone, and carbaryl, by directly or indirectly acting on baseline transmission of dopaminergic, serotonergic, and cholinergic receptors, respectively. Overall, differentiation in CINDA+ seems to produce higher numbers of these neuronal subtypes compared to differentiation in Electro medium, possibly due to the higher number of active electrodes these CINDA+ BrainSpheres produce in total. Haloperidol and buspirone bind to dopaminergic and serotonergic receptors, respectively. Haloperidol is an antagonist for the inhibitory dopamine D2 receptor [86,87]; thus, we observe an increasing effect of this drug on the wMFR of neuronal units. Interestingly, all NN contained more units responding with a decreased electrical activity, which is in line with the results of previous in vitro studies that observed only inhibitory reactions after acute treatment [126][127][128][129][130]. Görtz and colleagues suggested that this effect may occur due to a direct blockage of ion channels [126], which explains the rising number of units responding with decreased activity at the highest haloperidol concentration. Buspirone's primary MoA is binding to presynaptic inhibitory 5-HT 1A receptors as an agonist [88], thus producing an inhibitory action as we observe for most units within the NN in vitro. Carbaryl inhibits the enzyme acetylcholine esterase, thereby leading to an accumulation of choline in the synaptic cleft [89]. However, the second MoA of this insecticide is binding to nACh receptors [90,91,131]. Interestingly, previous in vitro studies only showed decreased electrical activity after exposure to carbaryl [132][133][134][135]. In this study, all NN exhibit neuronal units that respond in both directions, thereby covering both MoA. This might be due to an abundance of nACh receptor expression independent of the cholinergic synapse on GABAergic neurons [136].
Finally, we used the BrainSpheres in combination with spike sorting for setting up a test method, the human multi-neurotransmitter receptor (hMNR) assay for acute neurotoxicity testing that aims at enlightening the neurotoxic MoA for unknown test compounds in the future. As a small proof-of-concept study, we exemplified the use of this test method by studying the effects of the compounds TMT and emamectin with previously-described MoA for glutamatergic and GABAergic neurons. The hMNR confirmed the two MoA of the test compounds: (1) the enhanced glutamatergic activity by TMT-induced glutamate release [93]; (2) the reduced GABAergic neurotransmission caused by the GABA A and GABA C agonist emamectin [92]. Long-term exposure studies showed that TMT also affects synaptic vesicle fusion and recycling [137,138]. This could be a possible explanation for the wMFR reduction at higher TMT concentrations. Although the effect of emamectin was strongest in the GABAergic units, it was also observed in the unsorted and glutamatergic units. This is probably due to the high abundance of GABAergic units in these particular BrainSpheres, with 80 to 95% of the units reacting to GABA and bicuculline.
As this is a rather restricted case study, respective proof-of-concept studies for the applicability of the hMNR assay also have to be devised for the other neuronal subtypes. However, this small setup already demonstrates a high sensitivity for the two model compounds by detecting the compounds' effects in the nM range. In the end, we envision a test method setup that identifies all five different neurotransmitter receptors in the first identification phase followed by compound exposure with unknown substances. Spike sorting will identify compounds' MoA through this effort and hence deliver a neurotoxic MoA profile for each tested compound. The advantage of this system is that one analyzes neuronal units in a mixed neuronal/glia network context; however, information on the individual neuronal level is assessed.
Additional applications for acute substance testing in hiPSC-based BrainSpheres combined with spike sorting analyses are disease modeling and drug development. The pathophysiology of several neurological disorders such as Rett syndrome, autism spectrum disorders, schizophrenia, Down syndrome, and fragile X involves, amongst others, a disrupted GABA switch during brain development leading to an inhibitory/excitatory imbalance [139][140][141][142]. Moreover, they are suited as Parkinson's disease model and for untargeted disease modeling, revealing, to date, unstudied disease mechanisms and gene-environment interactions [62,143]. Another application can be envisioned in drug development for safety or efficacy evaluation, e.g., for seizure liability assessment [68,[144][145][146]. In addition, interference of compounds with neurotransmitter systems might also be an indication of their developmental neurotoxicity (DNT) potential. Therefore, in the future, this test method might also be a valuable addition to the current DNT in vitro testing battery [147] as test methods for substances' effects on neuronal subtypes were identified as one gap in current NAM-based DNT evaluation [148].

Summary and Conclusions
Taken together, we generated six different BrainSphere models by combining 2D neural induction protocols with 3D differentiation methods and showed distinguished neural differentiation patterns, although all protocols were based on dual SMAD inhibition. This emphasizes the importance of thorough characterization of each cell model and highlights the difficulties in comparing studies that use different media compositions. The different gene and protein expressions regarding neural subtype and receptor expression were also reflected in the functionality of the BrainSpheres measured on MEAs. To overcome the mixed signals of different neuronal subtypes, we applied spike sorting, which allowed us to distinguish between glutamatergic, GABAergic, dopaminergic, serotonergic, and cholinergic responses. Finally, this led us to introduce the hMNR assay which has possible applications beyond acute neurotoxicity for DNT testing, including in the field of disease modeling and for safety and efficacy evaluation in drug development.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cells12091270/s1, Figure S1: Immunocytochemistry staining of human induced neural progenitor cells (hiNPCs); Figure S2: Comparison of electrical activity of BrainSpheres (without 3D differentiation) for 7 weeks on microelectrode arrays (MEA); Figure S3: Comparison of electrical activity of 1 week 3D differentiated BrainSpheres for 7 weeks on MEA; Figure S4: Comparison of electrical activity of 2 weeks 3D differentiated BrainSpheres for 7 weeks on MEA; Figure S5: Spike raster plots (SRP) and microscopy images of 3 weeks 3D differentiated BrainSpheres differentiated for further 3 weeks on MEA; Figure S6: Modification of electrical activity by acute pharmacological modulation with glutamate and ap5/nbqx; Figure S7: Modification of electrical activity by acute pharmacological modulation with bicuculline; Figure S8: Modification of electrical activity by acute pharmacological modulation with picrotoxin; Figure S9: Modification of electrical activity by acute pharmacological modulation with haloperidol; Figure S10: Modification of electrical activity by acute pharmacological modulation with buspirone; Figure S11: Modification of electrical activity by acute pharmacological modulation with carbaryl; Figure S12: Specification of glutamatergic and GABAergic units for acute neurotoxicity testing; Table S1: Coatings used for neural differentiation of BrainSpheres generated with different neural induction protocols; Table  S2: Information about used compounds and their solvents; Table S3: Primer sequences; Table S4: p-values resulting from statistical analyses of the gene expression data.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.