Reverse engineering of feedforward cortical-Hippocampal microcircuits for modelling neural network function and dysfunction

Experimental design

Experimental timelines including a schematic of the microfluidic MEA design can be seen in Fig. 7. Our neural cultures comprised either commercially available rat embryonic cortical and hippocampal neurons (hereafter referred to as cortical-hippocampal networks), or lateral entorhinal cortex layer II neurons (LEC-LII), Dentate Gyrus granular neurons (DG-gl), CA3 pyramidal neurons (CA3-pyr) and CA1 pyramidal neurons (CA1-pyr), freshly dissected from adult McGill-R-Thy1-APP AD-model rats or APP/PS1 AD-model mice (hereafter referred to as adult entorhinal-hippocampal AD networks). The method for the dissection and culturing of adult hippocampal neurons was refined from our previously reported protocol for adult lateral-most LEC LII (lLEC-LII) neurons from APP/PS1 model mice43. For all cultures, coating and establishment of an astrocytic feeder layer was conducted using the same protocol.

Fig. 7
figure 7

Experimental timelines. (A) Four-nodal microfluidic chips were interfaced with either glass coverslips for immunocytochemistry (ICC) or microelectrode arrays (MEAs) for electrophysiology. All interfaces were coated with Poly-L-Ornithine and laminin, and a feeder layer of rat cortical astrocytes were plated two days prior to plating of commercially available rat cortical (node 1) and hippocampal (nodes 2-4) neurons. Spontaneously evoked activity was recorded every fourth day from 12 to 32 DIV, and electrical stimulations were induced after the recordings at 28 and 32 DIV. (B) For perturbed networks, all steps were identical up to 23DIV. At this point, amyloid beta (A(beta)) fragments were added to the cortical node. ICC was conducted 2, 5 and 9 days after perturbation to assess the change in concentration of amyloid beta oligomers in all nodes. Recordings were furthermore conducted 1, 5 and 9 days after the perturbation to evaluate the impact on the network functionality. (C) To reconstruct adult entorhinal-hippocampal circuits within the four-nodal microfluidic chips, brains of adult mice and rats were extracted and horizontally sectioned before dissection of the lateral entorhinal cortex layer II, and hippocampal subregional layers (dentate gyrus granular layer, and CA3- and CA1 pyramidal layers). Dissected entorhinal and hippocampal tissues were next dissociated into single cells enzymatically and manually by trituration. LEC LII, DG-gl, CA3-pyr and CA1-pyr were plated in nodes 1-4, respectively. A horizontally cut brain section viewed under a stereoscope using fiber optic illumination for contrast with delineation of entorhinal and hippocampal subregions is illustrated below the timeline. Lateral entorhinal cortex layer II (LEC LII) and hippocampal dentate gyrus granular layer (DG-gl), CA3 pyramidal layer (CA3-pyr) and CA1 pyramidal layer (CA1-pyr) are highlighted in blue.

Design & fabrication of microdevices

The four-nodal microdevices were designed using Clewin 4 (WieWeb Software, Enschede). Each node was 5 mm in diameter and ({60},{upmu }{hbox {m}}) high. Furthermore, the nodes were connected by 20 microtunnels. To induce unidirectional axonal outgrowth between the nodes, a geometrical design inspired by the Tesla valve was implemented within the microtunnel architecture24. All tunnels were ({350},{upmu }{hbox {m}})  long, ({10},{upmu }{hbox {m}})  wide and ({5},{upmu }{hbox {m}})  high. The design also had spine structures on the postsynaptic side to misguide axons trying to enter the microtunnel inlets. 50 nanoporous platinum electrodes of ({50},{upmu }{hbox {m}}) diameter were positioned evenly spread across the four nodes. Additionally, 9 electrodes were positioned within the microfluidic tunnels to confirm active connections between the nodes. Individual reference electrodes were embedded in the glass substrate within each of the four chambers on the platforms, with all four reference electrodes wired together and connected to the same contact pad on the MEAs. The design, as well as an image of a fabricated microdevice can be seen in Fig. S5. Fabrication of all microdevices was conducted according to our recently reported protocol24.

Coating of culturing platforms

Prior to coating, all microdevices were sterilized in UV light overnight in a biosafety cabinet. Consecutively, the distilled (DI) water in the devices was replaced by Poly-L-Ornithine solution (PLO) (Sigma-Aldrich, A-004-C) at a concentration of ({0.1},{hbox {mg/mL}})  and incubated at (37^{circ }hbox {C}), 5 % (hbox {CO}_2) for 2 h or overnight at ({4}^{circ }{hbox {C}}). Subsequently, all PLO was discarded and the microdevices rinsed three times for 10 min with milli-Q (MQ) water. After the last rinse, laminin solution consisting of ({16},{upmu }{hbox {mg/mL}})  natural mouse laminin (Gibco™, 23017015) in phosphate-buffered saline (PBS, Sigma-Aldrich, D8537) was added and the microdevices incubated at ({37}^{circ }{hbox {C}}), 5 % ({hbox {CO}_2}) for 2 h. A hydrostatic pressure gradient was established during all coating to ensure flow of the coating solution through the microtunnels. This was achieved by filling the chambers with an unequal amount of coating solution.

Seeding of astrocyte feeder layer

For the astrocytic feeder layer, a solution consisting of DMEM, low glucose (Gibco™, 11885084) supplemented with 15 % Fetal Bovine Serum (Sigma-Aldrich, F9665) and 2 % Penicillin-Streptomycin (Sigma-Aldrich, P4333) was prepared. Next, the coating solution was replaced by the astrocyte media, and rat cortical astrocytes (Gibco™, N7745100) were seeded at a density of ({100},{hbox {cells/mm}}^{2}), i.e., 2000 cells per microchamber. The astrocytes were allowed to expand for 48 h, before seeding of either embryonic cortical and hippocampal neurons or adult entorhinal and hippocampal neurons.

Embryonic neurons: plating and maintenance

Neural media consisting of Neurobasal Plus Medium (Gibco™, A3582801) supplemented with 2 % B27+ (Gibco™, A358201), 1 % GlutaMAX (Gibco™, 35050038) and 2 % Penicillin-Streptomycin (Sigma-Aldrich, P4333) was prepared. Rock Inhibitor (Y-27632 dihydrochloride, Y0503, Sigma-Aldrich) was additionally added to the media during the first two days of neural growth at a concentration of 0.1 % to increase neural survival. Rat cortical neurons from Sprague Dawley rats (Gibco, A36511) were plated at a density of ({1000},{hbox {cells/mm}}^{2}), equaling 20 000 cells in the first node of each microfluidic interface. The cortical neurons were thereafter allowed to settle in an incubator for 3 h before plating the hippocampal neurons. Next, rat hippocampal neurons from Sprague Dawley rats (Gibco, A36513) were plated at a density of ({1000},{hbox {cells/mm}}^{2}) in the three remaining nodes of each microfluidic interface. Half the neural media was replaced with fresh neural media 4 h after plating, and again after 24 h. From here on, half the neural media was replaced every second day throughout the experimental period.

Viral transductions

Viral transductions for inducing ubiquitous expression of fluorescent markers were used to verify unidirectional structural connectivity between the nodes, according to our previously reported method24. Viral vectors were prepared in-house at the Viral Vector Core Facility, NTNU. At 11 DIV, neural networks (n = 4) were transduced with AAV 2/1 serotype viruses loaded with either pAAV-CMV-beta Globin intron-EGFP-WPRE-PolyA or pAAV-CMV-beta Globin intron-mCherry-WPRE-PolyA under a CMV promoter to ubiquitously express GFP (nodes 1 and 3) or mCherry (nodes 2 and 4). 3/4 of the media in nodes 1 and 3 were removed to create a hydrostatic pressure gradient across the nodes, and viruses for expression of GFP were added at a concentration of (5e^2) viruses/cell (i.e. (1e^7) viruses per node). The cells were subsequently incubated for 3 h at ({37}^{circ }{hbox {C}}) and 5 % (text {CO}_2). Next, nodes 1 and 3 were filled all the way up with cell media, and 3/4 of the media in nodes 2 and 4 was removed before adding viruses to express mCherry at a concentration of (5e^2) viruses/cell. After 3 h of incubation, these nodes were also filled all the way up with media. Imaging of the networks were conducted at 28 DIV.

Amyloid beta perturbations

At 23 DIV, the cortical node in six four-nodal networks were perturbed using A(beta)-fragments. Due to the study’s complexity, experimental data were collected in two separate rounds: one for the control networks (n = 6) and another for the perturbed networks (n = 6). For both rounds, the same batch of astrocytes and neurons were utilized, but from different vials provided by the vendor (Gibco). 1 mg Beta-Amyloid Peptide (1-42, human) (ab120301, abcam) was dissolved in ({100},{upmu }{hbox {L}}) 1 % ({hbox {NH}_4{hbox {OH}}}) (458680025, Thermo Scientific). The solution was furthermore diluted to 2 mM in neuronal media. This concentration was chosen based on previous literature79, and initial concentration testing (Figs. S6 and S7), to induce A(beta)-oligomerization while still retaining network viability over prolonged periods of time. Subsequently, half the cell media in the cortical node was replaced by 2 mM A(beta)-solution. Immunocytochemistry was conducted 2, 5 and 9 days after the perturbations to assess the change in concentration of A(beta)-deposits in all nodes.

Animals models and genotyping

All animal experiments comply with the ARRIVE guidelines80, were approved by the Norwegian Food Safety Authority and carried out in accordance with the EU (European Union) Directive 2010/63/EU for animal experiments. The animals were held in standard lab cages (up to 5 animals per cage), with temperatures of 22 ± ({2}^{circ }{hbox {C}}), and kept at a light/dark cycle of 12:12 hours with access to food and water ad libitum. For optimization of the protocol we used a total of 43 adult McGill-R-Thy1-APP rats, either homozygote (+/+), heterozygote (+/-) or negative genotype (-/-), including 22 female and 21 male. Of these, we collected data from 14 rats (8 female and 6 male) and from a total of 2 female, adult APP/PS1 mice. See Supplementary Table S2 for a total overview of sex, age and genotype of all animals.

The APP/PS1 mouse model is a transgenic familial AD-model, with a C57BL/6J genetic background, co-expressing mutated human amyloid precursor protein (hAPP) (KM670/671NL) and mutated human presenilin 1 (L166). Both mutated proteins are driven by expression cassettes under the Thy1 promoter. This causes increased levels of A(beta) leading to formation of extracellular cortical A(beta)-deposits starting already from around 6 weeks of age81. The APP/PS1 mice were genotyped using a KAPA-kit as described in Hanssen et al.43. The McGill-R-Thy1-APP rat model is a transgenic familial AD model with a Wistar (HsdBrl:WH) genetic background, co-expressing mutated hAPP (APP K670M671delinsNL) and (APP V717F)82. The McGill-R-Thy1-APP rats were genotyped using quantitative PCR (qPCR) and genomic DNA isolated from ear tissue as described in previous studies83. Both homozygous and heterozygous rats were included in the study, in addition to controls.

Adult neurons: dissection

All tools were sterilized by autoclaving before use. Both mice and rats were deeply anesthetized using 5 % isoflurane gas (Abbott Lab, 05260-05) and checked for absence of reflexes before the head was decapitated with a guillotine for rats or surgical scissors (FST, 14007-14) for mice. Subsequently, dissection and sectioning of the brain were conducted as described previously for both mice and rats. See Hanssen et al. for a detailed description of the procedure43. Brains were extracted in a petri dish filled with Hanks Balanced Salt Solution (Thermo Fischer Scientific, 88284) kept on ice. Extracted brains were sectioned horizontally on a Leica VT1000 S vibratome (Leica Biosystems) (thickness of ({300},{upmu }{hbox {m}}), frequency of 5.3 Hz and speed of 8 mm/s) with 0.2 % AGAR (VWR, 20767.298) supporting the caudal end of the brain. Sectioned brain slices were immediately transferred to wells of a sterile 24-well plate held on ice. Sections were further transferred to a petri dish held on ice and viewed under a stereoscope with fiber optic lamps, allowing for sufficient contrast of the cytoarchitectonic landmarks. Dissection of the lateral entorhinal cortex layer II (LEC-LII) was conducted as described in our recent study for both mice and rats43. Additionally, the dentate gyrus granular layer (DG-gl), CA3 pyramidal (pyr) and CA1-pyr layers were dissected (Fig. 7C). The DG-gl can be separated from the molecular layer and hilus using its densely packed layer of granular neurons. Similarly, the pyramidal layer in CA3 and CA1 consists of densely packed pyramidal neurons that can be distinguished from superficial layers lucidum, radiatum and lacunosum-moleculare and deep layer oriens33. All dissected strips of tissue were placed by region in four separate sterile 15 mL tubes filled with dissection media (HABG-P) consisting of Hibernate-A (Thermo Fischer Scientific, A1247501) supplemented with 2 % B27+ (Gibco™, A3582801), 2.5 % GlutaMAX (Gibco™, 35050061) and 1 % Penicillin-streptomycin (Sigma-Aldrich, P4333), placed on ice.

Adult neurons: plating and maintenance

All dissected tissue was rinsed for debris three times by adding and removing 3 mL of HABG-P. Following this, all media was discarded and the tissue was dissociated for 15 min at ({37}^{circ }{hbox {C}}), 5 % ({hbox {CO}_2}), 20 % ({hbox {O}_2}) in Neural Isolation Enzyme Papain (Thermo Fischer Scientific, 88285) reconstituted in Hank’s Balanced Salt solution (Thermo Fischer Scientific, 88284) as described by the manufacturer (Thermo Fischer Scientific, MAN0011896). Subsequently, all dissociation enzyme was discarded, ({100},{upmu }{hbox {L}}) HABG-P was added and the tissue was manually triturated 10 times by using a ({100},{upmu }{hbox {L}}) pipette tip. Further, HABG-P was topped up to a total of ({3},{upmu }{hbox {L}}) and the vials of tissue were centrifuged at x200 g for 2 min in room temperature. After centrifugation, HABG-P was discarded and ({100},{upmu }{hbox {L}}) neural media was added before manual trituration 10 times using a ({100},{upmu }{hbox {L}}) pipette tip. The neural media consisted of Neurobasal Plus Medium (Gibco™, A3582801) supplemented with 2 % B27+ (Gibco™, A358201), 1 % GlutaMAX (Gibco™, 35050038), 2 % Penicillin-Streptomycin (Sigma-Aldrich, P4333), 0.1 % Rock Inhibitor (Y-27632 dihydrochloride, Y0503, Sigma-Aldrich) and 10 % Fetal Bovine Serum (Sigma Aldrich, 12106C). For adult neurons derived from mice, 0.1 % BDNF (Neurotrophins, 450-02) was added to the neural media. For adult neurons derived from rats, 0.1 % FGF2 (Thermo Fisher Scientific, 13256-029) was added to the neural media. Neurons were plated at a density of ({750},{hbox {cells/mm}}^{2}), equaling 15 000 cells in each node of the microfluidic chip. Half the neural media was replaced every second day throughout the experimental period.

Immunocytochemistry

For immunocytochemistry, both cortical-hippocampal and adult entorhinal-hippocampal neurons were plated in microfluidic platforms bonded to glass coverslips (VWR International, 24×24 mm No. 1 Menzel-Gläser). Fixation of cortical-hippocampal cultures was conducted using glyoxal solution based on the protocol by Richter et al.84. This fixative consisted of 20 % ethanol absolute (Kemetyl, 100 %), 8.0 % Glyoxal solution (Sigma-Aldrich, 128465) and 1 % acetic acid (Sigma-Aldrich, 1.00063) in MQ-water. Fixation of adult entorhinal-hippocampal cultures was conducted using PBS (Sigma-Aldrich, D8537) containing 0.4 M phosphate buffer and 4 % freshly depolymerized paraformaldehyde (PFA; Sigma-Aldrich, P6148). The cultures were washed with PBS to remove debris, before the fixative was applied for 15 min at room temperature. Subsequently, all chambers were washed with PBS 3 times for 15 min each. Next, 0.5 % Triton-X (Sigma-Aldrich, 1086431000) diluted in PBS was applied to permeabilize the cells. All chambers were again washed twice with PBS before a blocking solution consisting of 5 % goat serum (Abcam, ab7481) diluted in PBS was added to the cultures, and the cultures incubated at room temperature on a shaking table at 30 rpm for 1 h. Primary antibody solutions were prepared in PBS with 5 % goat serum, and antibody concentrations according to the ones listed in Table 1. Cultures were placed on a shaker table at 30 rpm at ({4}^{circ }{hbox {C}}) overnight. Following this, the primary antibody solution was removed, and the cell cultures were rinsed three times with PBS for 15 min each. Next, a secondary antibody solution consisting of 0.2 % secondaries and 5 % goat serum diluted in PBS was added to the cultures. The cultures were left to rest on a shaker table at 30 rpm for 3 h at room temperature. Prior to applying the secondary antibody solution, the prepared solutions were centrifuged at 6000 rpm for at least 15 min to remove precipitates. Subsequently, the secondary antibody solution was replaced by Hoechst (Abcam, ab228550) at a concentration of 0.1 % diluted in PBS, and the cultures left on a shaker table for another 30 min. Eventually, all cultures were washed three more times in PBS, then twice in MQ water prior to imaging.

Table 1 Antibodies and concentrations used for Immunocytochemistry.

All images were acquired using an EVOS M5000 microscope (Invitrogen). The microscope was equipped with DAPI (AMEP4650), CY5 (AMEP4656), GFP (AMEP4651) and TxRed (AMEP4655) LED light cubes and Olympus UPLSAP0 4X/0.16 NA and 20x/0.75 NA objectives. Post-processing of images was conducted in ImageJ/Fiji or Adobe Photoshop 2020.

Quantification of amyloid beta deposits

For quantification of amyloid beta deposits, ImageJ/Fiji was utilized. For each Day Post Perturbation (DPP; 2, 5, and 9), a total of four images were used for quantification in both perturbed and control networks, except for DPP 2 and 5 for control networks, where only two images were used. Each image was captured with a 20x objective, covering an area of 2046 pixels, equivalent to ({2803},{upmu }{hbox {m}}). All images were converted to 8-bit and thresholded using maximum entropy. Furthermore, the images were binarized, and particles separated using the watershed function. The built-in particle analysis function was subsequently used to count the number and size of the A(beta)-seeds, and only particles above ({5},{upmu }{hbox {m}}) were included in the final analysis. The threshold was set this high to assure only aggregates of amyloid beta were included in the analysis, and not noise caused by antibody precipitates or autofluorescence. This provided the best comparison of amyloid beta levels between networks stained at different time points.

Imaging and quantification of neurite length

For quantification of neurite length, phase contrast images acquired using a Zeiss Axio Vert V.1 brightfield 20x/53 WD/0.4 NA objective with an Axiocam 202 mono were used. For analysis, Neurolucida (MFB Bioscience) was used to quantify neurite length of the adult neurons.

Electrophysiological recordings

A MEA2100 workstation (Multichannel Systems) was utilized for all recordings with a sampling rate of 25000 Hz. An external temperature controller (TC01, Multichannel Systems) was used to maintain the temperature at ({37}^{circ }{hbox {C}}). All cultures were allowed to equilibrate at the workstation for 5 min prior to recordings, and recordings lasted for 15 min for cortical-hippocampal networks and 10 min for adult entorhinal-hippocampal neural networks. A 3D-printed plastic cover with a gas-permeable membrane was used to keep the neural networks in a sterile environment during the recordings.

Electrical stimulations

Stimulations of cortical-hippocampal neurons were performed following the recordings at 28 and 32 DIV using a train of 60 spikes at ± 800 mV (positive phase first) of ({200},{upmu }{hbox {s}}) duration with an interspike interval of 5 s. The most active electrode measured in spikes per second in the cortical cell population was chosen for stimulations. This was done to assure that the electrode chosen for stimulation had sufficient coupling with a neuron to induce activity.

Data analysis

All data preprocessing and analysis was conducted in Matlab R2021b, and graphs were plotted using the matlab function linspecer85, based on the color palette colorBrewer86. A 4th order Butterworth bandpass filter was used to remove noise below 300 Hz and above 3000 Hz, and a notch filter to remove noise at 50 Hz from the power supply mains. All filters were run using zero-phase digital filtering. For spike detection, the Precise Timing Spike Detection (PTSD) algorithm developed by Maccione et al.87 was utilized. A threshold of 9 times the standard deviation was chosen for the cortical-hippocampal cultures, with a maximum peak duration of 1 ms and a refractory period of 1.6 ms. For the adult entorhinal-hippocampal cultures, a standard deviation of 7.5 was chosen due to the lower signal-to-noise ratio of the activity in these cultures overall. The SpikeRasterPlot function developed by Kraus88 was adapted and used for creating raster plots.

Burst detection was conducted using the logISI algorithm developed by Pasquale et al.44. A minimum of four consecutive spikes was set as a hard threshold for a burst to be detected, and the maximum interspike interval to 100 ms. Network bursts were detected using the logIBEI approach89, with a minimum of 10 % of all active electrodes required to exhibit bursting activity for a network burst to be classified. To assess whether activity originating in the cortical node propagated to the hippocampal nodes, the sequence of electrode activation during network bursts was analyzed. For a network burst to be classified as having propagated through all four nodes, activity had to be detected sequentially in a feedforward direction, starting from the cortical node and progressing through hippocampal nodes 1, 2, and 3. Coherence index, a metric for assessing network synchrony, was calculated as the ratio of the standard deviation to the mean of the instantaneous firing rate90.

After binning the data into 50 ms time bins, functional connectivity was analyzed using cross-correlation. For each calculation, the maximum cross-correlation coefficient was computed between the binned spike trains of each node and up to 500 ms time-shifted copies of the activity from other nodes, with a maximum lag of 10 bins (500 ms). The cross-correlation was normalized so that the autocorrelations at zero lag were equal to 1, using the coeff specification in the Matlab function xcorr. Cross-correlation was selected as the correlation metric to capture the correlation between nodes that are farther apart, where the propagation of activity might take longer (e.g., between the cortical node and hippocampal node 3). The cross-correlation within nodes was evaluated by calculating the correlation between two and two electrodes within an individual node, and then assessing the average value of these correlations to get a single correlation value for each node. For the correlation between nodes, the spike trains from all electrodes in each node were combined, before the correlation between the combined spike trains between nodes was evaluated. This approach was adopted to minimize variability between nodes caused by differences in the number of electrodes detecting activity, given the limited number of electrodes per node. Since extracellular electrodes predominantly capture signals from the perisomatic region of neurons (within 50-({100},{upmu }{hbox {m}}) of the soma9), a single neuron is typically detected by only one electrode on these platforms. Furthermore, the 50 ms bin size used in the analysis ensures that activity from multiple neurons near an electrode is averaged within the same time bin during network bursts. This averaging reduces the potential impact of multiunit activity on internodal correlations and, consequently, diminishes the necessity for spike sorting in this context.

To visualize key network features, graphs were plotted using the electrodes as nodes, and cross-correlation as the edges (i.e. lines connecting the individual nodes). Edges with weaker connectivity than 0.1 were removed from the graph representations to highlight only the strongest, direct connections. Firing rate (Hz) was furthermore used to represent the node color and pagerank centrality the node size. The Louvain algorithm was used to perform community detection, delineating the nodes into distinct communities based on which nodes were more highly interconnected with each other than the rest of the network91. This analysis was used to determine the maximized modularity value of the networks. A high modularity value indicates that the network has dense internal connections within nodes and sparse connections between different nodes. The node edge color of the graphs were furthermore used to showcase which nodes belonged to the same community. The modularity value was used to assess how easily the algorithm managed to classify the nodes into distinct communities. All graph theoretical measures were calculated using the Brain Connectivity Toolbox developed by Rubinov & Sporns92.

To remove stimulation artifacts, stimulation data was run through the SALPA filter developed by Wagenaar et al.93. Additionally, 15 ms of the filtered data was blanked following each stimulation time point. Subsequently, the data was binned into 20 ms time bins, and peristimulus time histograms of the average response of the stimulations in each node plotted.

Statistical analysis

Generalized Linear Mixed-Effect Models (GLMMs) were used to assess differences in network parameters between control networks and networks perturbed with A(beta)-fragments. The GLMMs were generated using SPSS version 29.0.0.0. The network type, i.e. perturbed networks from DIV 24 to 32 (i.e. after perturbation) or control networks (i.e. all other data points for both perturbed and unperturbed networks) were used as a fixed effect. Considering that both control and perturbed networks received identical treatment until the point of perturbation, it was considered more appropriate to compare post-perturbation data points from the perturbed networks (represented by the pink bar in the GLMM) to the data points from the control networks throughout the experimental period and the perturbed networks before perturbation (represented by the blue bar in the GLMM). The network age was furthermore used as a random effect, and the network feature of interest as a target. To fit the data, a gamma distribution with a log link function was chosen based on the Akaike information criterion. Sequential Bonferroni adjustment was selected for multiple comparison. All other statistical analyses were conducted using a two-sided Wilcoxon rank sum test with the Matlab function ranksum.