Development of an objective index, neural activity score (NAS), reveals neural network ontogeny and treatment effects on microelectrode arrays

Microelectrode arrays (MEAs) are valuable tools for electrophysiological analysis, providing assessment of neural network health and development. Analysis can be complex, however, requiring intensive processing of large data sets consisting of many activity parameters, leading to information loss as studies subjectively report relatively few metrics in the interest of simplicity. In screening assays, many groups report simple overall activity (i.e. firing rate) but omit network connectivity changes (e.g. burst characteristics and synchrony) that may not be evident from basic parameters. Our goal was to develop an objective process to capture most of the valuable information gained from MEAs in neural development and toxicity studies. We implemented principal component analysis (PCA) to reduce the high dimensionality of MEA data. Upon analysis, we found the first principal component was strongly correlated to time, representing neural culture development; therefore, factor loadings were used to create a single index score—named neural activity score (NAS)—reflecting neural maturation. For validation, we applied NAS to studies analyzing various treatments. In all cases, NAS accurately recapitulated expected results, suggesting viability of NAS to measure network health and development. This approach may be adopted by other researchers using MEAs to analyze complicated treatment effects and multicellular interactions.

Principal component analysis of MEA parameters reveals temporal correlation, allowing for neural activity score derivation. Given the complexity and multivariate nature of the data, PCA was performed to reduce dimensionality and allow for easier visualization. After standard score normalization, all of the aforementioned parameters at all time points were included as data points for PCA. Examining the twodimensional projection of the first two principal components revealed a distinct pattern in the data (Fig. 2a). Adding a dimension of time (via colormap), this pattern was revealed to be a temporal separation of the data points, especially along the first principal component. Statistically, linear regression analysis supported this temporal component, as principal component 1 (PC1) is strongly correlated to time ( Fig. 2b; R 2 = 0.5441, p < 0.0001), indicating recapitulation of network ontogeny and maturation. After confirmation of this relationship, factor loading values for PC1 were examined to determine which factors (MEA parameters) contributed most strongly to this component. While substantial contributions were observed for many parameters, the strongest metrics were burst percentage, network burst percentage, number of spikes per burst, number of bursting electrodes, number of spikes per network burst, and synchrony index (Table 1). Notably, mean firing rate, the most common parameter analyzed in MEA studies, was the 11th-strongest contributor. Finally, these factor loading values were used to develop an individual index score-NAS (Eq. 1; see "Methods"). As NAS represents all aspects of neural network activity, it allows for assessment of neuronal network ontogeny and evaluation of the effects of various perturbations, such as stimulation, pharmacological treatment, or alternative culture conditions, which can typically be difficult to analyze if various parameters do not exhibit unidirectional changes. Additionally, NAS reduces the high variation often observed in individual MEA parameters, as evidenced by lower coefficient of variation for 24/25 (96%) measured parameters ( Supplementary Fig. S1).

Enhancement of neural network ontogeny is easily quantified via NAS. Regression analysis
served as initial support that NAS accurately measures neural network ontogeny, but we also sought to experimentally validate NAS in several conditions to further confirm this recapitulation. For initial validation, several experiments were performed to analyze enhanced neural network ontogeny and activity in response to different conditions known to enhance neural activity-namely, optimized culture media 13 and muscle-conditioned media treatment 14 . To examine the effects of optimized culture media, mixed neural cultures (HBG3-derived) were grown on MEAs in two different media conditions: DMEM/F12 & Neurobasal-based medium (DMNB) or BrainPhys-based medium (BP). While DMNB has traditionally been widely used to culture HBG3-derived and other neural cell lines, BP was developed for electrophysiological applications due to a more physiologically relevant formulation, resulting in increased electrophysiological function of various cell lines 13 . However, BP has not been evaluated on HBG3-derived neural cultures. In both DMNB and BP groups, the neurons began showing activity at approximately day 5, increasing over 3 weeks, as expected; however, the cells cultured in BP exhibited enhanced activity and network development, as indicated by the significantly higher NAS ( Fig. 3a; p < 0.0001, two-way repeated measures ANOVA).
To examine the effects of conditioned media on network ontogeny, mixed neural cells were treated with muscle cell (C2C12)-conditioned media (CM), which has previously been shown to significantly accelerate network activity and development 14 . Likewise, NAS analysis showed similar results and provided simple quantification ( Fig. 3b; p < 0.0001, two-way repeated measures ANOVA) of this accelerated network development.
Disruption of neural network ontogeny is easily quantified via NAS. In addition to measuring neural network activity enhancement, we also sought to validate NAS on more complex culture conditions and for quantifying disruption of network activity. Microglia, the resident immune cells of the central nervous system (CNS), are being increasingly implicated in neurodegenerative diseases and have been shown to be neurotoxic in many conditions [15][16][17][18] ; therefore, we decided to explore co-culturing microglia with mixed neural cultures on To examine whether this disruption is contact-dependent or the result of secreted factors, neural cultures were treated with BV2-conditioned media at 10 days (similarly to the co-culture experiment described above). Similar to the co-culture condition, BV2-conditioned media treatment also disrupted network function (Fig. 4b), suggesting a role for microglia-secreted factors in neural network disruption. To examine whether this effect was exacerbated by microglial activation, BV2 cells were stimulated with two concentrations of the pro-inflammatory endotoxin lipopolysaccharide (LPS; 10 ng/mL and 100 ng/mL) for 24 h prior to conditioned media collection. LPS-stimulated BV2-conditioned media disrupted network function in a concentration-dependent manner, with unstimulated BV2-conditioned media causing significant disruption (p < 0.0193, two-way mixed model, Tukey's post-hoc test), followed by 10 ng/mL LPS stimulation (p = 0.0003) and 100 ng/mL LPS stimulation (p < 0.0001), providing further support for NAS as a viable method to quantify complex treatment effects and evaluate disruption of electrophysiological function.
NAS summarizes neural activity for neurotoxicology screening. Advances in MEA technology 8,19 have led to adoption of MEAs for neurotoxicological screening 20 Given the potential of NAS to consolidate many functional MEA parameters, we sought to determine its applicability to neurotoxicity screening.
For this analysis, NAS values were calculated from MEA toxicity screening of 52 compounds from the NTP or ToxCast libraries 21,22 (Fig. 5a-c). In previous studies 21,22 , the authors performed a network formation assay (NFA) using primary cortical neurons, measuring 17 parameters of activity in response to compound treatment over 12 days on MEAs to determine compound effects on network formation. Additionally, viability testing was performed to measure cytotoxicity. For each of these assays, EC 50  Of the 52 compounds we analyzed, 33 were found to have measurable effects in the network formation assay (defined as a decrease in activity by 3× median absolute deviation from control) for at least one activity parameter (though the specific parameter(s) differed among compounds), and 26 compounds were found to

Discussion
Advances in MEA technology, including multi-well MEA plates, incubated recording setups, and constantly improving software, allow for higher throughput than previously possible 18,19 , though analysis has traditionally been limited to simple parameters, primarily mean firing rate. Only recently have researchers begun incorporating advanced metrics of network activity in these screening approaches [20][21][22] . These advanced metrics have provided researchers with tools to record from entire neuronal populations and analyze complex neuronal network dynamics. Multi-well MEAs enable high-throughput neuronal recordings, facilitating their adoption for drug/ toxicity screening applications and evaluation of complex culture conditions. However, the information that can be gleaned from MEAs has been hampered by limited analytical methods and tools, as well as high variation. Here, we present a novel method to overcome these limitations-development of an index, neural activity score, that incorporates and consolidates traditional MEA measurements into a single quantitative value that can be used to objectively evaluate neuronal network development and function across various culture conditions, treatments, and neural cell sources. This is valuable not only for basic neuroscience research on neuronal networks, but also translational research and preclinical studies. To facilitate adoption, further validation, and improvement, codes used to calculate NAS can be found at the following open source repository: https:// doi. org/ 10. 5281/ zenodo. 39393 10. Despite these advantages, NAS is not free from its own limitations-particularly those associated with its derivation from principal component analysis. As NAS is essentially a linear weighted combination of individual activity parameters, parameters changing in equal (weighted) magnitude but in opposite directions could negate each other, resulting in no net change. While this example demonstrates the potential for different electrophysiological profiles with identical NAS, the overall goal of NAS is to provide a holistic interpretation of neural development and maturity, suggesting that these "equal but opposite" changes still represent cultures with similar overall states of maturation. For example, a culture with increased firing rate but decreased synchrony may represent an overall similar level of maturity-especially when comparing different cell types or culture conditions that may exhibit variations in ontogeny. NAS provides a way to "balance" these changes that would be difficult Another potential limitation is sensitivity to detect meaningful changes compared to noise, which is difficult to evaluate. Technologically, this sensitivity is determined by the MEA platform, recording settings, and analysis thresholds, such as the thresholds for determining spikes and bursts. For NAS, this sensitivity relies on the particular changes in activity. In one case, minute changes in multiple parameters in the same direction are summed, which would be measured as a larger overall change in NAS than the individual changes; however, in cases such as the one described above, small changes in opposite directions may go undetected. While this case may result in relatively low sensitivity, it also reduces variation, providing additional confidence that detected changes represent meaningful development and maturation.
One final limitation involves incorporation into existing analytical pipelines. While we encourage use and modification of the provided code and formulas to calculate NAS, we recognize that these are additional elements and steps to add to already complex pipelines and software packages that vary across individual labs and equipment. These extra steps may hinder rapid adoption for users with different MEA platforms and custom analyses, especially in other programming languages (e.g., MATLAB, R).
The results presented here demonstrate the value of NAS to assess potential developmental neurotoxicity (DNT) hazards, a field with a widely recognized need for more sensitive, less variable, and higher throughput functional assays [21][22][23][24] . The mixed neural cultures used for NAS derivation and the primary cultures analyzed in the network formation assay are both maturing networks, derived from embryonic stem cells or isolated from neonatal rodents, respectively. As a result, NAS is well-suited for analysis of maturing neural networks, as is necessary in DNT studies, covering a range from non-active to full maturity, with synchronized network bursting. The application to developing networks from multiple cell sources suggest NAS has substantial value for improving the use of MEAs for toxicity screening and drug development.
As the concern over drug development costs continues to rise, scientists are noticing several recurring problems, including the reproducibility crisis and inadequacies of current screening assays, in vivo models, and other preclinical studies [25][26][27] . For neural assays, specifically, assays have traditionally used simple endpoints such as viability and morphological analysis (i.e. neurite outgrowth) for screening, primarily due to scalability 28 . However, electrophysiological endpoints are often more sensitive and allow for assessment of electrophysiological toxicity, which involves separate-and highly time-sensitive-mechanisms 8,29,30 . By improving result www.nature.com/scientificreports/ interpretation, NAS will facilitate incorporation of functional measures into screening programs focused on cytotoxicity and morphology. Index scoring has been used extensively in clinical settings and in vivo; for example, neurological deficits in amyotrophic lateral sclerosis (ALS) and Parkinson's disease (PD) are assessed using the Revised ALS Functional Rating Scale (ALSFRS-R) 31 and United Parkinson's Disease Rating Scale (UPDRS) 32 , respectively. Stroke severity is frequently measured using various scales [i.e., modified Rankin scale (mRS) 33,34 and NIH stroke scale (NIHSS)] 35 , and these have been shown to correlate strongly with patient outcomes and be useful for therapeutic evaluation 34 . The simplified analytical pipeline provided by these indices is vital to detecting effects (or lack thereof) in clinical and preclinical studies. Due to this, a need has been recognized to develop multivariate approaches and index scores for in vitro approaches, as well 36,37 . Similar analysis pipelines provided by index scores could be especially valuable for screening assays, allowing for improved hit detection when screening potential neurotoxicants or therapeutics. Several composite scores have been developed to condense information from multiple toxicity assays for specific compound classes (e.g. endocrine disruptors, halogenated aliphatics), previously 38,39 . Here, we developed NAS using a similar approach to condense the high-dimensional data from MEA recordings into a single measurement with reduced variation that can be used to easily and consistently evaluate compound effects on neural activity, as opposed to analyzing many different parameters individually. www.nature.com/scientificreports/ This reduced variation and improved interpretation could help identify and/or narrow down compounds to examine and develop further, saving time and money wasted on poor candidate compounds. Likewise, improved in vitro studies could help reduce the necessity of in vivo studies, which are expensive, time-consuming, and have ethical and practical concerns due to a myriad of potential endpoint measurements and species differences that can contribute to high variability and difficulty determining true treatment effects 40 .
The validation studies presented here indicate that the NAS formula provides an easily interpretable measure of neural network health/functionality and overall effects of perturbation. By compiling all MEA metrics as opposed to individual metrics (i.e., mean firing rate), NAS represents all aspects of neural network function, which can provide more consistent analysis and results interpretation/reporting. Additionally, NAS has the potential to provide increased sensitivity over a collection of individual parameters, as NAS was more sensitive than the individual parameter average for 61.5% of compounds. This result was interesting, demonstrating the utility of implementing relative parameter weights. Since NAS was derived from how all parameters contribute to development/maturation, this result indicates that this approach may describe treatment effects in a more holistic manner than analysis of individual parameters alone, which only describe certain aspects of activity. However, when specific parameters are of interest, we suggest incorporating NAS as an additional metric for screening, not as a complete replacement, as a summary statistic for electrophysiological function and neural network maturation. Additionally, larger training data sets and/or other optimization may allow for improved sensitivity in the future.
Two of the three compounds for which NAS was found to be most sensitive, MPP+ and picoxystrobin, share similar toxic mechanisms, both inhibiting mitochondrial electron transport chain complexes 41,42 . While further research would be needed to determine if this is more than a coincidence, it does suggest mitochondrial function as a sensitive predictor of neurodegeneration. This finding supports a wealth of evidence linking mitochondrial dysfunction to neurodegenerative diseases, in some cases prior to symptom onset and diagnosis [43][44][45] . Using NAS to analyze and compare various compound classes in more detail may allow for deeper insight into toxicity mechanisms for different compound classes or varying therapeutic potential in drug discovery studies. To this point, the compounds analyzed here were primarily organohalogens, halogenated alkanes, and pyrethroids. Future analysis of additional compound classes with varying properties and mechanisms of actions (e.g., organophosphates) would provide further validation for NAS, as well as additional data on the sensitivity of MEAs and electrophysiological screening to assess developmental neurotoxicity. Indeed, a larger screen of 27 organophosphates detected neurotoxic effects for approximately 1/3rd of tested compounds 21 ; however, further studies are needed to make comparisons to the classes examined here and whether other assays may provide higher sensitivity for specific classes of compounds.
Lastly, challenges in analyzing complex and large data sets have been widely acknowledged across multiple assays and techniques, including high-throughput screening, image analysis, and flow cytometry [46][47][48][49][50] . These challenges include high variability, difficulty interpreting results across multiple metrics, and reproducibility-problems that are only exacerbated when examining complex/emergent phenomena that may be difficult to quantify otherwise, such as neuronal network function. While we developed and validated NAS using MEA data, many of the solutions posed for the aforementioned techniques also utilized PCA and other dimensionality reduction methods, suggesting a similar index scoring approach may be useful for these, and other, applications to gain a deeper understanding of important results.
BV2 microglia cells (gift from Dr. Jae-Kyung Lee, University of Georgia, Athens, GA) were cultured according to previously published protocols 51 . Briefly, cells were thawed and seeded on tissue culture-treated plates at approximately 5-10,000 cells/cm 2 and passaged at 60-80% confluency. Cells were maintained with media changes every other day with neural medium consisting of DMEM/F12 (ThermoFisher) supplemented with 5% fetal bovine serum (FBS) (GE Healthcare, Chicago, IL), 2 mM l-glutamine (ThermoFisher), and 1% penicillin/  www.nature.com/scientificreports/ streptomycin (ThermoFisher). For lipopolysaccharide (LPS) treatment, cells were treated with 10 ng/mL or 100 ng/mL LPS in neural medium for 24 h before conditioned media was collected and centrifuged to remove any cells or cellular debris.
MEA preparation, recording, and data processing. 48-well MEA plates (Axion Biosystems) were prepared according to manufacturer's protocol. Briefly, plates were coated with 0.1% polyethyleneimine (PEI) (Sigma) for 1 h at 37 °C, rinsed with sterile water, and allowed to air dry in a biosafety cabinet overnight. The following day, plates were coated with 20 µg/mL laminin (Sigma) for 2 h at 37 °C prior to cell seeding. Mouse neural cultures (see above) were seeded and allowed to adhere for 1 h, then maintained in full neural culture media (see above) supplemented with GDNF (Peprotech) and CNTF (Peprotech) (10 ng/mL each) with media changes every 3-4 days throughout the 3-week recording period. Neuronal activity was recorded using the Maestro system (Axion Biosystems) and AxIS software v2.1-2.5 (Axion Biosystems) with the following settings: band-pass filter (Butterworth, 300-5000 Hz), spike detector (adaptive threshold crossing, 8xSD of RMS noise), burst detector (100 ms maximum inter-spike interval, 5 spikes minimum, 10 spikes minimum for network bursts, 10 ms mean firing rate detection window). Recordings were performed daily for 2 min at 37 °C after allowing plates to acclimate to the Maestro system.
Raw data files were processed offline using the Statistics Compiler function in AxIS. Statistics Compiler output files were processed in Microsoft Excel (Microsoft Corporation, Redmond, WA) and with custom Python scripts to organize and extract individual parameter data for each well of each MEA plate and for data normalization.
Raster plots generated with Neural Metric Tool v1.2.3 software (Axion Biosystems): https:// www. axion biosy stems. com/ produ cts/ softw are/ neural-module# appli catio ns Neural activity score calculation. After initial data processing, normalization (to z-score values), and outlier removal (− 3 > z > 3), JMP 14 (SAS Institute, Cary, NC) was used to conduct principal component analysis. All parameters (Supplementary Table S1) were used for all wells (replicates) at 5-19 days in vitro (DIV). The first two principal components were used to visually analyze the temporal separation of the data (Fig. 2a), then the first principal component was used for linear regression analysis to determine the extent of correlation to time. Finally, the factor loadings for the first principal component were calculated to show the extent of contribution for each individual MEA activity parameter (Table 1).
Factor loadings for principal component 1 were then implemented as coefficients in a formula incorporating, and ultimately consolidating, all of the measured individual MEA parameters into an individual index score-NAS-defined as the sum of each measured parameter value multiplied by its factor loading value for each well (replicate) at each time point (Eq. 1).
where β i are the factor loading values and x i are the z-normalized measured values for each parameter.
Analysis of DNT hazard screening data. Raw MEA data [*.raw files generated via the Maestro system and AxIS software (Axion Biosystems), see above] from previous studies 21,22 was processed through the same analysis pipeline described above. Additional processing for neurotoxicity data was based on methods described by Shafer et al. 21 , including area under curve (AUC) calculation, Hill function fitting, and EC 50 extrapolation. Specifically, AUC values for each compound and concentration were calculated in Python 3 using the trapezoidal rule (numpy.trapz() function) to integrate normalized NAS values over time (see Data Availability below for more information about custom Python codes). Concentration-response curves (NAS AUC vs. concentration) were generated via nonlinear least squares regression ([Inhibitor] vs. normalized response model) in GraphPad Prism 8.2.0 (GraphPad Software Inc., San Diego, CA) for each compound with Hill slope = − 1.0, and EC 50 values were extrapolated from the resulting curves.
EC 50 values corresponding to cytotoxicity (Supplementary Table S2) that were used for sensitivity analysis were reported from previous studies 21,22 . EC 50 values for NAS and average MEA parameter were calculated as described above.