Modulating electrophysiology of motor neural networks via optogenetic stimulation during neurogenesis and synaptogenesis

Control of electrical activity in neural circuits through network training is a grand challenge for biomedicine and engineering applications. Past efforts have not considered evoking long-term changes in firing patterns of in-vitro networks by introducing training regimens with respect to stages of neural development. Here, we used Channelrhodopsin-2 (ChR2) transfected mouse embryonic stem cell (mESC) derived motor neurons to explore short and long-term programming of neural networks by using optical stimulation implemented during neurogenesis and synaptogenesis. Not only did we see a subsequent increase of neurite extensions and synaptophysin clustering, but by using electrophysiological recording with micro electrode arrays (MEA) we also observed changes in signal frequency spectra, increase of network synchrony, coordinated firing of actions potentials, and enhanced evoked response to stimulation during network formation. Our results demonstrate that optogenetic stimulation during neural differentiation can result in permanent changes that extended to the genetic expression of neurons as demonstrated by RNA Sequencing. To our knowledge, this is the first time that a correlation between training regimens during neurogenesis and synaptogenesis and the resulting plastic responses has been shown in-vitro and traced back to changes in gene expression. This work demonstrates new approaches for training of neural circuits whose electrical activity can be modulated and enhanced, which could lead to improvements in neurodegenerative disease research and engineering of in-vitro multi-cellular living systems.

such as potentiation and homeostatic plasticity, have made it difficult to control neural circuit development for in vitro applications 8 . Furthermore, studies have shown varying degrees of emerging phenomena that are correlated to neuronal plastic events that would imply "learning" or "memory storage" in neural circuits [9][10][11][12][13] These varying degrees extend from neural mechanisms that regulate potentiation to feedback mechanisms with respect to associations to stimulation. Currently, in-vitro micro electrode array (MEA) systems have aimed to study the effects of stimulation protocols in the resulting network dynamics for the past 30 years [13][14][15][16][17][18][19] . While some of these studies have induced transient changes in the recorded signals, most have focused on changes of number of action potentials in mature primary neurons and in learning as a result of feedback to strengthen associations to the stimulation patterns. Due to the fact that these phenomena emerge from multiple pathways related to neuronal plasticity, finding new approaches to enhance and modulate these plastic responses in a long-term way would be highly advantageous in the field.
Because most of these studies on plasticity have focused on modulation of mature neurons, we hypothesized that we could induce plasticity-related long-term electrophysiological changes in in-vitro neural networks, by implementing training regimens during early stages of differentiation, i.e. neurogenesis, coupled with training regimens during network formation, i.e. synaptogenesis. To this end, we used 465 nm pulsed light to induce depolarization in specific temporal patterns to implement a training regimen on differentiating channelrhodopsin (ChR2)-expressing MEBs grown in suspension and continued the training regimen after seeding them on functionalized glass to allow for neurite extension and network development. In parallel, we seeded the MEBs on MEA chips while continuing training regimens to be able to characterize the network dynamics during network development. We first characterized responses to the optogenetic stimulation during neurogenesis by assessing morphological parameters, followed by analyzing the resulting electrophysiological responses including the network's synchronicity, firing patterns in power spectra and the system's responsivity to stimulation during recording. These results showed a direct correlation between perturbations during differentiation and plastic responses occurring during network formation. Finally, through RNA sequencing studies, we observed genetic changes that serve to explain the observed modulations in these neuronal systems.

Materials and methods
Experimental model and subject details. mESC culture and differentiation. A feeder layer of mitomycin-C inactivated mouse embryonic fibroblast (MEF) was seeded at a cell density of ~ 3.5 × 10 4 cells/cm 2 and cultured in DMEM (Dulbecco's modified Eagle's medium) supplemented with 1% fetal bovine serum, 1% Lglutamine and 1% penicillin-streptavidin. Subsequently, HB9:GFP transgenic mouse embryonic stem (mES) cells transfected with (Channelrhodopsin) ChR2-TdTomato were seeded on the feeder layer at a ratio of 1.5 mES cells per MEF. Media was changed to mESC proliferation medium and replaced daily.
Differentiation begun prior to mESC colonies coming in contact with one another. Cultures were trypsinized (0.05% Trypsin) after being exposed to embryonic stem cell differentiation medium (eDM) for an hour, and later seeded on 100 mL low adhesion dishes in 10 mL eDM. The next day, floating cells were collected to separate the culture from adhered non-neuronal lineages. On the following day, embryoid bodies (EBs) were replated in eDM supplemented with 1 µM of retinoic acid (RA) (Sigma Aldrich, MA) and 1 µM puromorphine (PM) (STEMCELL Technologies, MA). On D5, EBs were resuspended in fully supplemented eDM (FS eDM) supplemented with RA and PM plus 10 ng/mL of growth factors, glial cell-derived neurotrophic factor (GDNF) and ciliary neurotrophic factor (CNTF) respectively. Media was changed daily.
MEA fabrication. Platinum micro-electrode array chip were fabricated on borofloat glass wafers following standard lithographic techniques. Photoresist LOR3A (MicroChem, MA, USA) was spun on clean substrates, followed by spin coated layer of photoactive S1805 (Dow Chemical Comp., MI). Spun substrates were exposed using an EVG 620 (i-line) aligner (EV Group Inc., Tempe, AZ) and developed. Ti/Pt (1:3) was later evaporated for a total thickness of 1,000 Angstroms. To passivate traces between the detection area and contact pads, 300 nm of silicon nitride was deposited on the entire substrate using Plasma-Enhanced Chemical Vapor Deposition System (PlasmaLab International, WA). Finished chips were diced into 49 × 49 mm squares and fitted with acrylic wells bonded by Dow Corning Sylgard 184 (Ellsworth Adhesives, IL).
To ensure that thermally generated noise voltages were below membrane voltage fluctuations during neuronal electrical activity, electrode impedance was decreased through electrochemical treatment of the electrodes 20 . Platinum black was formed on top of clean electrodes to achieve high surface area using a Gamry Reference 600 Potentiostat. The galvanostatic deposition was achieved by running a chronopotentiometry experiment at 2.83E-6 A/cm 2 versus Ag/AgCl for 15 s in a solution of dihydrogen hexachloroplatinate (0.08 mM H 2 PtCl 6 -6H 2 O, Sigma Aldrich, with 0.25 g/L of (CH 3 COO) 2 Pb Alfa Aesar) for a total of 21.45 ng (1.71E-2 ng/µm 2 ) of crystallized platinum. The impedance reduction due to the platinum black deposition was examined by electrochemical impedance spectroscopy.
Electrophysiology recording. MEA measurements were performed using a MEA 2,100-Lite Amplifier (Multi Channel Systems MCS GmbH, Germany) at 37 °C. Electrical activity from ChR2 + MEBs cultured on MEA were measured every other day. Measurements were performed in dark at a sampling rate of 10 kHz for 20 min in FS eDM with sealed covers to keep CO 2 concentration stable.
Optogenetic stimulation was performed using a laser diode (LD) Driver (Doric Lenses, Quebec, Canada) attached to a single LD, Blue 465 nm with a Fiberoptic Patchcord for an incident intensity of 10 mW/mm 2 which ensured that the intensity at the samples was still above the 1 mW/mm 2 limit for ChR2 activation even after refraction from the lid and media. Because embryoid bodies were below 200 µm of diameter, the incident intensity of the blue light could penetrate across the entire EB 21 Potentiation regimens consisted of 5 ms pulses at  20 Hz for one second every other second for 1 h for a total of 1,800 cycles, while stimulation during recording  consisted of the same pattern for only 10 cycles at equidistant times during the 20 min recording sessions. Measurements and calculations of heat transfer from beam to the media showed that the energy imparted would cause temperature fluctuations of less than 1 °C, which are less than common fluctuations from normal cell handling.
Immunostaining and imaging. Immunocytochemistry was performed on samples fixed with 4% paraformaldehyde and treated with 0.05% Triton-X. Permeabilized samples with blocked with 4% bovine serum albumin at 4 °C overnight. Samples were later stained with primary antibodies at 4 °C overnight, followed by staining of secondary antibodies at room temperature for 2 h. All antibodies were diluted in Antibody Diluent (ThermoFisher Scientific, MA). The used antibodies were 1) anti-NeuN antibody for neural populations, (clone A60, Alexa Fluor 555 conjugate), 2) anti-GAD65/67 polyclonal goat antibody (Santa Cruz Biotechnology), 3) anti-vGlut polyclonal guinea pig antibody (Synaptic Systems, Germany), and 4) DAPI as a nuclear stain. For synaptophysin clustering quantification between MEBs, the used primary antibody was anti-synaptophysin 1 polyclonal chicken antibody (Synaptic Systems, Germany). It was counter stained with goat anti-chicken IgY AlexaFluor 647 (Abcam, MA), respectively, along with DAPI as a nuclear stain. After washing overnight at 4 °C in blocking buffer, samples were mounted with Prolong Diamond antifade (ThermoFisher, MA) and imaged using Zeiss 880 Confocal microscope (Carl Zeiss Microscopy).
Scanning electron micrographs. SEM images were taken at 1 kV after grounding the entire array through the contact pads, to avoid charging the insulating Nitride layer, at 100-300X magnification.
RNA extraction. Samples of MEBs were carefully collected and centrifuged at 14,000 rpm for 15 min. After aspirating supernatant, samples were flash frozen using liquid nitrogen and immediately stored at -80 °C, until mRNA extraction. Total RNA was collected and purified using the RNeasy Mini Kit Part 1 (Qiagen). The total RNA concentration was quantified using a Nanodrop spectrophotometer.
Quantification and statistical analysis. Spike/burst analysis. Multi-Channel Analyzer software was used for counting and analytically extracting temporal parameters of fast events. Raw data was digitally filtered using a 2nd Order Butterworth high pass filter (cutoff frequency: 200 Hz). Action Potentials (APs), were detected as "spikes" by setting a threshold at 5 × standard deviations from the noise magnitude distribution. Analysis of firing rate behaviors and burst parameters followed spike detection 22 . For this experiment, burst detection was defined by the following parameters: • Max interval to start burst: 50 ms • Max interval to end burst: 50 ms • Min. interval between bursts: 100 ms • Min. duration of burst: 50 ms • Min. number of spikes in burst: 4 This analysis was done for active electrodes which were defined as electrodes which recorded at least 10 AP/ min. Spectral analysis. Spectral analysis was performed on the slow component of raw data to assess modulations in network behavior. Data extracted as ASCII files were filtered in MATLAB using Butterworth 2 nd order high pass digital filter with a cutoff frequency of 200 Hz. This filtered data was then used to detect the occurrences of spikes. Spectral components in the frequency domain from this binned spike data were obtained through Fast Fourier transform (FFT) between 0.1 and 200 Hz to remove DC components from data and detect frequency components occurring in bursts or clusters of spikes.
FFT was obtained for subsequent non-overlapping intervals of 10 s across the initial 4 min of spontaneous activity and normalized by the area under the curve. Calculated spectra were smoothed using a 3-point window moving average. Data was stored in 2D matrices, summed across electrodes and averaged across MEAs.
Synaptophysin cluster counting. Stack images were superimposed for 10 microns in ImageJ. Then a binary threshold was set so that only saturated pixels were conserved. Images were sampled ten times at regions in between MEBs with a 150-by-150 µm area across three biological repeats by day.
Quantifying network synchronicity. Overall network correlation was assessed through the automated use of spike data and a customized MATLAB code that calculated cross-correlation for discrete functions, as follows: Scientific RepoRtS | (2020) 10:12460 | https://doi.org/10.1038/s41598-020-68988-y www.nature.com/scientificreports/ Here, f and g denote the discreet functions at point n, which represent the recorded MEA signals. Furthermore, f* denotes the complex conjugate and m is the lag or displacement, meaning a feature in f occurring at n that occurs in g at n + m 23 . Spike train data for each channel was cross-correlated with every other channel. Results were normalized so the autocorrelations have unit value at zero lag. The value at zero lag (t ch-x = t ch-y ) was stored for each correlation in a 60-by-60 matrix relating every channel to each other and plotted as a heat map and a bar graph was used to plot the average.

RNA sequencing analysis.
Each sample that was categorized as an N for the RNA Sequencing consisted of an entire culture of mESC grown and differentiated from a 35-mm well which resulted in around 100-150 MEBs. Biological triplicates were obtained by thawing and evenly dividing three different vials of 2E6 mESC (< P7) evenly across three wells for a total of 9 seeded wells, and each was handled independently. After expansion and neural induction steps (at D2), one well pertaining to each of the 3 thawed vials were snap frozen. Furthermore, one well pertaining to each of the 3 thawed vials were assigned to be trained, and the remaining three served as controls. At D9, all MEBs were collected in their separate samples and snap frozen for RNA extraction and downstream RNA Sequencing.
The RNA-seq libraries were prepared with Illumina's 'TruSeq Stranded mRNA-seq Sample Prep kit' (Illumina). The libraries were quantitated by qPCR and sequenced on one lane on a HiSeq 4,000 for 101 cycles from one end of the fragments using a HiSeq 4,000 sequencing kit version 1. Fastq files were generated and demultiplexed with the bcl2fastq v2.20 Conversion Software (Illumina). Adapters and low-quality bases were trimmed from reads using Trimmomatic 24 (v0.36) with parameters LEADING:28 TRAILING:28 MINLEN:30. The trimmed reads were quasi-mapped to Gencode's M19 transcriptome using Salmon 25 (v 0.8.2) with additional parameters -seqBias -gcBias -numBootstraps = 30. Transcript expression value were summarized to the gene-level and corrected for average transcript length using the "lengthScaledTPM" method 26 . TMM-normalized 27 log2 counts per million (cpm) values (prior.count = 3) were calculated and only genes with > log2(0.5 cpm) in at least 3 samples were analyzed using the limma-trend method 28 . Three pairwise comparisons were made (C vs D2, S vs D2 and S vs C) and multiple testing correction was done separately for each comparison using the False Discovery Rate method 29 ; significant differential expression was indicated at FDR p-value < 0.05.
Statistical analysis. All data sets from electrophysiological responses were extracted from raw recordings using Matlab (Mathworks, Natick MA, USA). To validate the MEB-derived neural networks by measuring the response to known signaling molecules, one MEA per signaling molecule was used, selecting 15 active electrodes. Testing responses of neurite extension to training regimens during neurogenesis was done by monitoring and imaging 20 separate embryoid bodies seeded on gridded coverslip and measuring and averaging at the 4 cardinal points for each measured embryoid body. Response of presynaptic puncta to stimulation regimens was measured for 10 images from 40x objectives of separate cultures. Electrophysiological responses were measured for three cultures grown on separate MEAs per experimental group, with each replicate value obtained from averaging active electrodes from each culture. All results, unless stated otherwise, are expressed as mean + /-SEM. Statistical comparisons were performed in OriginLab using one-way ANOVA and repeated measures ANOVA (when appropriate) for group and longitudinal difference test and followed by post-hoc Tukey test for multiple comparisons when applicable. Differences were considered statistically significant if p < 0.05.

Results
Optogenetic stimulation was used on mESC-derived MEBs to implement training regimens during two important stages of neural development: neurogenesis (while still in suspension) and synaptogenesis (seeded on functionalized glass or MEAs) (Fig. 1a). Training regimens consisted of periodic stimulation with 5 ms pulses at 20 Hz in 1 s intervals for an hour ( Supplementary Fig. 1a). This regimen has been shown to enhance axonal growth 30 , and thus would suggest that it could lead to a shift in structural potentiation in a neural network. The regimen was repeated every 24 h as differentiation occurred within the EBs, with an expectation that consistent repetition would enhance the potentiation and cause long-term changes in the firing patterns of the network. Following established differentiation protocols of mESC towards mature motor neurons 31-33 , the described training regimen was started at D2 of differentiation, at which point stem cells have been induced towards neuronal lineages, and specialization and maturation of motor neurons has been shown to take place in the subsequent 7 days (Fig. 1b). Since one of the transcription factors that drove differentiation, retinoic acid, is light sensitive, media was changed every single day immediately after stimulation to ensure that stimulation effects on MEBs were not artifacts (i.e. false positives) caused by photodegradation of factors ( Supplementary Fig. 1b) 34 . Furthermore, since the differentiation was monitored with the expression of the motor neuronal marker Hb9 through a GFP reporter, we used the plateau of GFP expression between D8 and D9, as an indicator that D9 was an appropriate time point for seeding the MEBs on glass ( Supplementary Fig. 1c). Thus, after these 7 days (D2-D9) of differentiation, stimulated (S) and non-stimulated (NS) cultures were seeded on MEA chips (Fig. 1c). Careful seeding practices were applied to ensure that ~ 20 MEBs were seeded within the sensing area of the MEAs for a ~ 50% coverage by the MEBs (Supplementary Fig. 2). Seeding in this manner ensured empty space between clusters for the extension of processes, even though some nearby clusters would start fusing into larger clusters. The resulting two groups of samples seeded on MEAs were further subdivided into two more experimental groups, referring to whether or not a training regimen was continued during network formation on chip for the consequent 15 days (D10-D25). For ease of discussion, S or NS prior to a colon (e.g. S:X or NS:X) will refer to the presence or lack thereof of stimulation, during neurogenesis, while S or NS written after a colon (e.g. X:S or X:NS), indicates the presence or absence of stimulation during synaptogenesis (Fig. 1a).
Scientific RepoRtS | (2020) 10:12460 | https://doi.org/10.1038/s41598-020-68988-y www.nature.com/scientificreports/ The electrical activity of the resulting neuronal cultures was measured with the MEA system and the raw data was filtered to remove low frequencies (< 200 Hz), to remove undesired voltage artifacts (e.g. stimulation artifacts), and extract action potentials recorded as spiking events (Fig. 1d). A two-step procedure was used to remove false positives from the analyzed data: (1) the detection threshold was set at a value at which no positives would be detected from the ground electrode, then (2) the recorded spikes at each electrode were inspected to ensure that the detected spikes had the appropriate voltage phases relating to action potentials: depolarization, repolarization and refractory period.

MEB cultures form active neural networks with excitatory and inhibitory populations. In this
work, neural networks were cultured from intact MEBs, in contrast to growing them as a monolayer after dissociation. The long-term goal of our study is the modulation of electrical activity of the MEBs towards downstream implantation in in-vivo or in-vitro experimental systems and modulating the functionality of such systems through the resulting interaction. When cultured in their intact form, MEBs tend to keep their spheroid shape, while extending processes which contain neurites that form networks as they undergo synaptogenesis (Fig. 2a). Furthermore, dense web-like neurite structures form within the spheroid itself (Fig. 2b) and both excitatory (vGlut) and inhibitory (GAD65/67) receptors stain positively (Fig. 2c).
Network formation was validated by exposing MEB cultures grown on MEAs (Fig. 2d) to varying concentrations of commonly used exciting and inhibiting signaling molecules for 5 min: L-glutamate, acetylcholine, cyclic    (Fig. 3e). When monitoring the difference of pre-synaptic clusters per unit area at D13 and D15, the statistically significant difference indicated that optogenetic stimulation during neurogenesis evoked physiological responses on two important aspects of neural network development: neurite extension and presynaptic clustering (Fig. 3e).
MEB network synchronicity is amplified by stimulation during neurogenesis and synaptogenesis. Network synchrony is a common parameter used to characterize a developing neural network, as it gives information on the network's plasticity and connectivity. Various studies have successfully shown that the presence of chronic stimulation results in improved network synchrony [36][37][38] . In our study, we wanted to observe the long-term effects of stimulation regimens on the network synchrony and determine if these effects were amplified or shifted when the training regimen during neurogenesis was extended during synaptogenesis. From the raster plots of the spontaneous activity recorded at D21, the increased level of synchronous activity was notable between NS:S and S:S samples versus S:NS and NS:NS (Fig. 4a). This can be appreciated by the peaks above the raster plots, which correspond to a summation of the activity across all electrodes, where synchronous networks would result in discrete peaks whereas in samples that lacked coordinated firing, the resulting line plot seemed to lack any peaks. Similarity between electrode recordings was quantified with cross-correlation in order to quantify synchronous behavior. Values for the similarity across the network were obtained by calculating cross-correlation for all electrode combinations (Supplementary Fig. 3). For this analysis, only spontaneous recordings of active  (Fig. 4b). This showed that synchronous behavior extended across the entire network and was markedly higher for networks that were stimulated during synaptogenesis.
Interestingly, when the network wide mean synchronicity was calculated for each recording day, a trend of higher synchrony was observed for samples that had been exposed to some form of training regimen (NS:S, S:NS or S:S) but no statistical significance was observed at D11 (ANOVA, n = 3; Spectral density elucidates changes in steady state firing. Conventionally, electrophysiological behavior is characterized by firing rate during set epochs and burst parameters ( Supplementary Fig. 4). However, when analyzing these parameters during spontaneous firing, there was no discernable trend in the change of long-term firing rate or burst parameters between experimental groups. However, when observing the spike data during steady state of a more mature neural network (D25), there were deviations on how the spike firing clustered into bursts, despite the fact that no clear change in the number of spikes was observed (Fig. 5a). We accredited this seeming conflict between the quantitative and qualitative data to the selection method of the burst detection parameters (See Quantification and statistical analysis). In order to avoid arbitrariness in the selection of these parameters, we decided to characterize the data in the frequency domain. For this reason, we focused on characterizing spontaneous firing recorded on MEAs by comparing changes in the power spectrums of recorded signals calculated through Fourier transforms (Fig. 5b). To obtain spectral profiles, binned spike counts were divided into 10-s-long contiguous windows and transformed to the frequency domain, thus representing the power spectrum as a function of time (Fig. 5b). When initially calculating the power spectral density (PSD) and observing between the DC frequency and the Nyquist frequency, we noticed that most of the components appeared below 7 Hz for all samples. For this reason, we compared samples between 0.1 Hz (to remove DC component) and  (Fig. 5c). Moreover, the post hoc Tukey test indicated a statistically significant difference between power spectra values between NS:S Neurogenetic stimulation changes the opto-response of MEB networks. Another aspect of consideration on the effect of training MEBs during neurogenesis was whether the early stage perturbation had some effects on how the later-stage network would respond to the same perturbation. To study this, we recorded responses to optogenetic stimulation from sets of samples that had not undergone the training regimen during neurogenesis (Fig. 6a) and compared them to those set that had undergone such regimen (Fig. 6b). Initial observation showed a difference between how the networks responded when stimulated early in the network development (D11) versus more mature networks (D25). For example, when early networks, which had a low spontaneous firing rate (D11) were stimulated, there would be a very notable evoked response during stimulation followed by a quiescent state, where the network would barely fire before returning to the baseline spontaneous firing rate. In contrast, more mature networks (D25), would still show an evoked response during stimulation but would automatically return to baseline firing rate right after stimulation ceased. What was interesting was that the quiescent time after stimulation for early S:S networks were notably shorter than those from the NS:S samples (Fig. 6ab). Moreover, at D25, while NS:S samples would return to the same baseline firing rate right after stimulation stopped, S:S samples showed a transient change in firing rate for several seconds after the stimulation stopped ( Fig. 6a-b). www.nature.com/scientificreports/ To quantify this behavior, the evoked firing rate during stimulation (FR stim ) and the post-response firing rate (FR post ) were compared to the firing rate prior to stimulation (FR pre ) for the three instances of stimulation within recording for each of the three MEA networks for both experimental groups (Fig. 6c). While the fold-change increase of firing rate FR pre to FR stim decreased with time for both NS:S (repeated measures ANOVA with Greenhouse-Geisser correction, n = 3; F(1.48, 11.83) = 14.79, p = 1.12E-3) and S:S (repeated measures ANOVA with Greenhouse-Geisser correction, n = 3; F(1.88, 15.02) = 11.02, p = 1.31E-3 (because more mature networks would have a higher baseline firing rate), when comparing the amount of evoked action potentials during stimulation (FR stim /FR pre ), S:S samples seemed to respond more strongly to stimulation than NS:S samples (Fig. 6d). One-way ANOVA determined a statistically significant difference between NS:S and S:S FR stim /FR pre values for D13 (n = 9; Additionally, the quiescent state response post-stimulation observed in early days (D11, D13 and D15), reflected itself in FR post being less than FR pre , resulting in FR post /FR pre < 1 for NS:S and S:S samples. We observed that this transient decrease in firing rate was statistically significantly shorter for the S:S samples than the NS:S for D11 (ANOVA, n = 9; F(1,16) = 19.95, p = 3.9E-4) and D13 (ANOVA, n = 9; F(1,16) = 9.49, p = 7.2E-3) (Fig. 6e). Repeated measured ANOVA indicated that FR post /FR pre ratios increased for both NS:S (Greenhouse-Geisser corrected, n = 9; F(3.06, 24.48) = 36.92, p = 2.69E-9) and S:S (n = 9; F(7,56) = 5.66, p = 5.63E-5). Furthermore, at later days of network development, it was notable that FR post /FR pre was ~ 1 for NS:S, meaning that the steady state firing rate was indistinguishable from that immediately following the termination of stimulation. On the other hand, S:S samples showed FR post /FR pre values above 1 from D17 forward, indicating that the network would transiently increase in firing rate right after stimulation. One-way ANOVA showed that this increase between FR post /FR pre values for S:S and NS:S was statistically significant for D17 (n = 9; F(1,16) = 12.19, p = 3E-3), D21 (n = 9; F(1,16) = 6.94, p = 0.018) and D23 (n = 9; F(1,16) = 9.91, p = 6.23E-3), while only showing a non-statistically significant trend for D19 (n = 9; F(1,16) = 2.16, p = 0.16) and D25 (n = 9; F(1,16) = 3.76, p = 0.071). It is relevant to mention that these effects were observed while there was no perceivable change in efficiency of the blue light to activate the ChR2 ion channels and evoke a response in the networks (Supplementary Fig. 6). These observations were corroborated by repeated measures ANOVA performed at p < 0.05, which showed no statistically significance change in efficiency (repeated measures ANOVA, n = 12; F(2,22) = 1.25, p = 0.31).
To further study how the training regimens affected network response, we also quantified the evoked response reflected in the network's synchronicity for the initial stimulation done on the initial spontaneous interval of recording. For this purpose, raster-plots of the average values of cross-correlation (as calculated for the analysis in Fig. 4) were calculated using 10 s bins across the entire 20 min of recording (Fig. 6f). When quantifying the short term effect of stimulation during recording had on network synchronicity, by comparing  (Fig. 6g).
Changes evoked by stimulation during neurogenesis result in genetic changes. Given the effects on neurite extension, presynaptic clustering, frequency profiles and network response to stimulation that were observed as a result of the presence of training regimens on MEBs during neurogenesis, we proceeded to determine genetic changes that could provide possible mechanistic explanations. Total messenger RNA sequencing was performed and analyzed for stimulated (S) and non-stimulated (NS) MEBs at D9, as well as EBs at D2. The differentially expressed genes in MEBs that underwent training regimens during neurogenesis were compared to those that did not, both with respect to the genetic expression of EBs sampled prior to differentiation (at D2). A total of 749 differentially expressed genes between S and NS with p < 0.05 were detected and clustered and color coded with respect to the differential expression of D2 (Fig. 7a). There were 200 genes that were upregulated during control differentiation, but this upregulation was lessened for samples that underwent training regimen (black bar), while the upregulation of 172 genes was amplified for those same samples (red bar). On the other hand, there were 202 genes whose downregulation was stagnated for samples with training regimen (yellow bar). For 173 genes, the control downregulation was further amplified after stimulation (blue bar). Something important to note was that this observed differential expression did not include changes in phenotype populations, matching the immunostaining observations ( Supplementary Fig. 7). This indicated that training regimen during differentiation did not seem to noticeably disrupt the rate of phenotype specification or generation of the neural populations that generally result from the differentiation protocol (Table 1). This suggests that training regimens affected other functional pathways rather than altering the differentiation of populations. For further analysis, a more stringent threshold (p < 0.0005) was set to detect the most promising genes as key factors for the behavioral changes seen in stimulated MEB cultures. This threshold resulted in 97 differentially expressed genes for the black cluster (Fig. 7b), 63 differentially expressed genes for the red cluster (Fig. 7c), 77 differentially expressed genes for the yellow cluster (Fig. 7d) and 71 differentially expressed genes for the blue cluster (Fig. 7e). From this pool, a thorough literature study was used to identify gene targets that had been reported to be related to known neural development and function (Table 2, Supplementary Fig. 9).

Discussion
Engineering of neural circuits is a critical step for biomedicine as well as various neuroengineering efforts. This work aimed to expand upon previous findings on achieving long term modulation of a neural network's behavior. 10,13,16,[39][40][41][42] . Here, we showed a direct relation between changes in network function and development and training regimens implemented during neurogenesis and synaptogenesis. To our knowledge there has been no electrophysiological characterization of mESC-derived neural networks' electrical response where a relationship between training during neurogenesis and synaptogenesis is examined. Furthermore these effects were correlated and supported by RNA-seq studies, where we showed that the training regimen affected targets related to network development, rather than affecting the differentiation of neural populations. In mESC-derived embryoid bodies, stimulation during neurogenesis has shown that the system enters a heightened state for neurite extension, which could be correlated to enhanced network formation. Furthermore, our study suggests that the presence of stimulation during neurogenesis has caused the system to be more sensitive and responsive to external stimulation Figure 7. RNA Sequencing shows differential expression as a result of optical stimulation during neurogenesis. a. Heat map of standard deviation of differential expression for genes with p < 0.05 (n = 749). Genes were primarily clustered for: (1) genes that would overexpress during differentiation and underexpressed due to stimulation, (2) genes that would overexpress during control differentiation and overexpressed further due to stimulation, (3) genes that would underexpress during control differentiation and stimulation minimized that underexpression and (4) genes that would underexpress during control differentiation and stimulation amplified that underexpression. (first color column in order: black, red, yellow, blue). Significantly differentially regulated genes, with p < 0.0005 (n = 307) were extracted as column plots for: b. black, c. red, d. yellow and e. blue clusters.
Scientific RepoRtS | (2020) 10:12460 | https://doi.org/10.1038/s41598-020-68988-y www.nature.com/scientificreports/ after network formation, which could have implications on research studies where MEBs are integrated in in-vivo nerve tissues that have suffered injury. By implanting "trained" MEBs that have an enhanced responsiveness to stimulation, the therapeutic effects that have been observed could in turn be enhanced as well. 3 Furthermore, spectral analysis proved to be a useful tool in order to quantify changes in network behavior that might be missed by simply comparing firing rates. While firing rate did not clearly reflect any alteration to the steady state, transforming the electrical signal to the frequency domain served as an informative descriptor of alterations in the network dynamics as a result of the training regimens. While spectral analysis has already been proven to be effective in characterizing neural signals 43,44 , the results of this study further elucidates its ability to detect modulations in the activity of neural networks that might not come across clearly from temporalbased parameters. Using FFT to characterize the network's firing patterns as a response to the training regimens during network formation complements various other network effects that have been observed in this study. For example, it has been shown that optogenetic stimulation is capable of causing a higher level of synchrony in cultures consisting of primarily excitatory neurons 37 . This would suggest a more structured firing in the form of a more consistent frequency profile across time, since synchronous behavior would be more viable in structured firing patterns versus a sporadically firing network. It is relevant to point out that the frequency profile in the S:S cultures while not significant on its own, confirms a modulation in the electrophysiological activity of this in-vitro model in response to the presence of the regimen across the development of the neural network. Furthermore, it serves as a proof of concept for the modulation of electrophysiological activity using this type of experimental setup. Following the results shown here, it would merit to explore further alterations to the network behavior using other analogous training regimens across the different developmental stages and using both spectral and temporal parameters to characterize what changes arise.
Another very notable observation regarding how the presence of training regimen implemented during neurogenesis resulted in long-term changes in the neural network was the difference of how the neural system responded to the stimulation during recording. During early culture on chip (D11-D15), stimulation would cause the system to enter a quiescent state following the forced state of rapid firing from the stimulation, yet this quiescent state was statistically significantly shortened for samples that had been stimulated during neurogenesis. Moreover, these distinctions extended beyond early stages of network formation, which also manifested as prolonged transient increases of network synchrony. This is relevant because, while tetanic stimulation has been www.nature.com/scientificreports/ correlated to potentiation, no change in baseline firing rate was observed. Therefore, if no long-term changes had been caused by training regimen during neurogenesis, there should have been no reason for the optogenetically evoked firing rate on chip to have shown any differences to the NS:S counterparts. This enhancement of response in samples that were subjected to training regimens during neurogenesis points to long-term systematic changes on how the system would respond to the stimulation during downstream stages of network development.
To further validate these observations, RNA sequencing was performed on MEBs showing that a large number of genes were found to have been differentially expressed as a result of the training regimen. To focus our findings, the search was narrowed to genes showing a fold change with p < 0.0005 and key genes related to neural function were selected. Upon observation of this narrowed list, we found that following training regimen during neurogenesis the expression of key developmental genes Npdc1 and Crabp2 more closely resembles that of more mature tissue. Other differentially expressed genes implied an improved development of key neural aspects such as axonal and dendritic growth (Tuba1a) and nucleation of microtubules (Tubgcp4). Furthermore, stimulated MEBs seemed to overexpress important factors related to vesicular signaling such as Snap47, Vamp2, Reep2, and more importantly Cacng7. Furthermore, various genes are directly related to improved synaptic function and plasticity. For example, upregulation of Aplp1, Cnih2 and Insyn1 would be targets of interest that could account for the structured firing patterns and the changes in the network responses to evoked stimulation. Considering the change in expression of these genes related to vesicle transport, synaptic transmission and neural development, the observations of shorter quiescent times post-stimulation for samples that underwent training regimens could suggest that the training regimen enhanced the network for a quicker recovery from synaptic fatigue, which is an activity-dependent temporary inability of neurons to transmit signals 45 . Moreover, it is worth noting that Cnih2 and Insyn1 are regulators of inhibitory signaling related to the slowing down of signal transmission to improve its reliability, yet a lot is unknown about how these targets affect larger scale pathways and neuronal function. While many factors responsible for regulating network plasticity have been discovered 11,36,[46][47][48] , how these interconnect and result in a steady state firing pattern, and furthermore, how could this pattern be controllably shifted is still to be determined. By continuing plasticity studies which incorporate training with respect to developmental stages and coupling the changes observed in the electrophysiology behavior to changes in genetic expression, these questions might be better addressed. Moreover, expanding this study to obtain protocol-specific responses in gene expression would be crucial in understanding neural plasticity and learning mechanisms.
Nevertheless, the results presented in this study add a new dimension to neural circuit engineering, by taking advantage of the synergy between responses to stimulation during differentiation and the plasticity that emerges www.nature.com/scientificreports/ during network formation, and correlation to gene expression. This solidifies our initial hypothesis that implementing training regimens during early stages of neural circuit development should be considered as a critical factor for studying neural plastic phenomena. In order to further this study, improved training regimens and spatial stimulation patterns during neurogenesis and synaptogenesis could be used to further the characterization of the network response. Furthermore, this approach could be expanded to in-vivo systems in order to examine if these alterations in fact result in behavioral changes in a developing organism, and if the enhancement of mESCderived MEBs have an effect regeneration studies by enhancement recovery of nerve tissue injury.