Astrocytes, the most abundant non-neuronal cell type in the mammalian brain, are crucial circuit components that respond to and modulate neuronal activity via calcium (Ca²⁺) signaling. Astrocyte Ca²⁺ activity is highly heterogeneous and occurs across multiple spatiotemporal scales—from fast, subcellular activity to slow, synchronized activity across connected astrocyte networks—to influence many processes. However, the inputs that drive astrocyte network dynamics remain unclear. Here we used ex vivo and in vivo two-photon astrocyte imaging while mimicking neuronal neurotransmitter inputs at multiple spatiotemporal scales. We find that brief, subcellular inputs of GABA and glutamate lead to widespread, long-lasting astrocyte Ca²⁺ responses beyond an individual stimulated cell. Further, we find that a key subset of Ca²⁺ activity—propagative activity—differentiates astrocyte network responses to these two major neurotransmitters, and may influence responses to future inputs. Together, our results demonstrate that local, transient neurotransmitter inputs are encoded by broad cortical astrocyte networks over a minutes-long time-course, contributing to accumulating evidence that significant astrocyte-neuron communication occurs across slow, network-level spatiotemporal scales. These findings will enable future studies investigating the link between specific astrocyte Ca²⁺ activity and specific functional outputs, which could build a consistent framework for astrocytic modulation of neuronal activity.

Animals

Experiments were carried out using young adult mice, in accordance with protocols approved by the University of California, San Francisco Institutional Animal Care and Use Committee (IACUC). Animals were housed in a 12:12 light-dark cycle with food and water provided ad libitum. Animal housing rooms were kept at 68–74 degrees Fahrenheit and 30–70% humidity. Male and female mice were used whenever available. Transgenic mice used in this study were Cx43^fl/fl mice⁵⁶ from the Bhattacharya Lab (UCSF, USA) and EAAT2-tdT mice⁵⁷ from the Yang Lab (Tufts University, USA). For in vivo imaging, all experiments were performed at the same time each day. 

Surgical procedures

For viral expression for ex vivo experiments, neonatal Swiss Webster or C57Bl/6 (P0–3) mice were anesthetized on ice for 3 min before injecting viral vectors (AAV5.GfaABC₁D.GCaMP6f.SV40 [Addgene, 52925-AAV5], AAV9.hGfap.pinkFlamindo, pENN.AAV9.Gfap.iGluSnFr.WPRE.SV40 [Addgene, 98930-AAV9], or AAV5.GFAP(0.7).RFP.T2A.iCre [Vector Biolabs, 1133]). Pups were placed on a digital stereotax and coordinates were zeroed at lambda. Four injection sites in a 2 × 2 grid pattern over V1 were chosen. Injection sites were 0.8–0.9 mm and 1.6–1.8 mm lateral, and 0 and 0.8–0.9 mm rostral. At each injection site, 30–120 nl of virus were injected at a rate of 3 nl/s at two depths (0.1 mm and 0.2 mm ventral/below pia) using a microsyringe pump (UMP-3, World Precision Instruments).

For viral expression for the in vivo experiments, adult C57BL/6 mice (2–4 months at the time of surgery) were administered dexamethasone (5mg/kg, s.c.) >1 hour before surgery, and anesthetized using 1.5% isoflurane (Patterson Veterinary Supply, 78908115). After hair removal and three alternating swabs of 70% ethanol (Thermo Fisher Scientific, 04-355-720) and Betadine (Thermo Fisher Scientific, NC9850318), a custom-made titanium headplate was attached to the skull using cyanoacrylate glue and C&B Metabond (Parkell, S380). A 3mm craniotomy was made over the right visual cortex. Virus was injected at two sites in right visual cortex at coordinates centered on +2.4mm and +2.7mm medial/lateral, +0.35mm and +0.65mm anterior/posterior and -0.3mm dorsal/ventral from lambda. 300nL of AAV5.GfaABC₁D.GCaMP6f.SV40 (Addgene, 52925-AAV5) was injected at each site through a glass pipette and microsyringe pump (UMP-3, World Precision Instruments). After allowing at least ten minutes for viral diffusion, the pipette was slowly withdrawn and a glass cranial window implanted using a standard protocol.

Ex vivo two-photon (2P) imaging and uncaging

Coronal, acute V1 slices (400-µm thick) from P28–32 (bath-application) and P27–42 (uncaging) mice were cut with a vibratome (VT 1200, Leica) in ice-cold slicing solution containing (in mM) 27 NaHCO₃, 1.5 NaH₂PO₄, 222 sucrose, 2.6 KCl, 2 MgSO₄, 2 CaCl₂. Slices were transferred to pre-heated, continuously aerated (95% O₂/5% CO₂) standard artificial cerebrospinal fluid (ACSF) containing (in mM) 123 NaCl, 26 NaHCO₃, 1 NaH₂PO₄, 10 dextrose, 3 KCl, 2 MgSO₄, 2 CaCl₂. Younger mice were sliced in the same solutions for GCaMP bath application of LY379268 and Baclofen (P20–25), Pink Flamindo (P20–22), and GluSnFR (P14–17). Slices were kept at room temperature until imaging. Bath-application experiments were performed at room temperature and 2P uncaging experiments were performed at 29°C using an in-line heater (TC-324B and SH-27B, Warner Instruments). To block neuronal action potentials during all slice imaging experiments, except for GluSnFr recordings, TTX (1 µM) was added to the ACSF > 10 min before imaging and remained in the circulating bath for the duration of the experiments.

Images were acquired on an upright microscope (Bruker Ultima IV) equipped with two Ti:Sa lasers (MaiTai, SpectraPhysics). Laser beam intensities were modulated using two independent Pockels cells (Conoptics) and images were acquired by scanning with linear galvanometers. Images were acquired with a 16×, 0.8 N.A. (Nikon) or a 40×, 0.8 N.A. (Nikon) water-immersion objective via photomultiplier tubes (Hamamatsu) using PrairieView (Bruker) software. For GCaMP imaging, 980 nm excitation and a 515/30 emission filter were used. For RFP imaging, 980 nm excitation and a 605/15 emission filter were used. For Pink Flamindo and Alexa Fluor 594 imaging, 1040 nm excitation and a 605/15 emission filter were used. Images were acquired at 1.42 Hz frame rate, 512 × 512 pixels and 0.64–1.61 µm/px resolution. For GluSnFR imaging only, images were acquired at 6.21 Hz frame rate, 200 × 200 pixels and 0.64 µm/px resolution, with 980 nm excitation and a 515/30 emission filter.

For bath-application experiments, a 5-min baseline was recorded to monitor spontaneous activity, after which receptor agonists were added along with a fluorescent dye (Alexa Fluor 594 Hydrazide) to assess the time at which drugs reached the imaging field/recording chamber (except for Pink Flamindo due to spectral overlap). An ACSF washout period (> 10 min), followed by a TTX incubation period (>10min), occurred between trials when imaging the same slice sequentially for bath-application of different receptor agonists or uncaging of different RuBi-subtypes. To account for any changes resulting from prior agonist exposure or uncaging, we alternated the order of agonists between concentrations or RuBi-subtypes between slices.

For simultaneous 2P imaging and uncaging, a second Ti-Sa laser beam was tuned to 800 nm and controlled using an independent set of linear galvanometers from those used for scanning. Laser beam intensity was modulated using an independent Pockels cell (Conoptics) to achieve a power measurement of ~2–8 mW at the slice. The beam paths for imaging and uncaging were combined after the linear galvanometers using an 855-longpass dichroic mirror (T855lpxr, Chroma). The uncaging laser was calibrated each experimental day by burning spots into a fluorescent slide. RuBi- compounds (300 µM) and TTX (1 µM) were added to the ACSF >10 minutes before imaging each slice. Fields of view (FOV) were chosen based on the location of GCaMP expression, which was often biased to/brighter in deeper cortical layers (distance of FOV from pia: 615 ± 196µm [mean ± SD, n = 121 FOV]). Prior to imaging at each FOV, a 60-s period was recorded to identify potential uncaging sites. Areas of GCaMP expression that exhibited moderate levels of spontaneous Ca²⁺ activity were chosen as uncaging sites. For FOVs with sequential GABA/glutamate uncaging, a continuous 5-min recording was used to monitor activity in each FOV. For FOVs with three sequential rounds of glutamate uncaging, a continuous 12.5-min recording was used to monitor activity in each FOV. Each recording began with a 2.5-min baseline period and at the 2.5-min mark, neurotransmitter was uncaged with 10 x 100 ms pulses, 100 ms apart. Sequential recordings of GABA/glutamate uncaging within the same FOV were separated by > 20min. Rounds of sequential glutamate uncaging were separated by ≥ 25min. Voltage from the uncaging laser Pockels cell was recorded to mark the time of uncaging pulses. Because RuBi-GABA and RuBi-glutamate are light-sensitive, care was taken to ensure experiments were carried out in minimal light. The computer screen and red-shifted headlamp were covered with two layers of red filter paper (Roscolux #27 filter, Rosco) and all indicator lights on equipment were covered. 

In vivo 2P imaging

In vivo 2P imaging was performed on the same microscope as ex vivo imaging, via a Nikon 16x, 0.8 N.A. water-dipping objective with a 2x-optical zoom (frame rate: 1.7Hz, FOV: 412µm², resolution: 512x512 pixels). Animals were given > 1 week after surgery for recovery and viral expression. They were then habituated on a custom-made circular running wheel over at least two days, and for a cumulative time of at least 2.5 hours, before recording. After habituation, mice were head-fixed on the wheel and movements were recorded by monitoring deflections of colored tabs on the edge of the wheel using an optoswitch (Newark, HOA1877-003). To compute wheel speed, a detected break in the optoswitch circuit was determined when the absolute value of the derivative of the raw voltage trace was at least 2 standard deviations above the mean. For recordings with little movement (std < 0.1), this threshold generated false positives, so a set threshold of 0.1 was used. The number of breaks in the optoswitch circuit per second was then calculated, and using the circumference and number of evenly spaced colored tabs at the edge of the wheel, the wheel speed was determined and used for all subsequent analyses using speed. Movement periods were defined by wheel speed ≥ 10 cm/s and movement bouts that were separated by ≤ 2 s were considered one event. To ensure that movement related dynamics were not included in stationary analysis, data was excluded from < 10 s around identified movement periods. GCaMP was imaged with 950nm excitation light and a 515/30 emission filter. Recordings lasted 30 minutes.

Ex vivo pharmacology

The following concentrations of each pharmacological reagent were used for experiments as indicated in the text: Tetrodotoxin Citrate (TTX, 1 µM, Hello Bio); Carbenoxelone disodium (CBX, 50 µM, Tocris Bioscience); R(+)-Baclofen hydrochloride (5–100 µM, Sigma-Aldrich); (1S,3R)-ACPD (t-ACPD, 5–100 µM, Tocris); LY 379268 disodium salt (100 µM, Tocris); Alexa Fluor 594 Hydrazide (0.1–2 µM, ThermoFisher Scientific); RuBi GABA trimethylphosphine (RuBi-GABA-Pme₃, 300 µM, Tocris); RuBi-Glutamate (300 µM, Tocris); CGP 55845 hydrocholoride (10 µM, Tocris); and LY 341495 (10 µM, Tocris). 

Immunohistochemistry and image quantification

After recording, slices from 2P imaging experiments were immersed in 4% PFA for 30 min and switched to 30% sucrose for one day at 4°C before being embedded in OCT and stored at -80°C. Slices were re-sectioned coronally at 40 µm on a cryostat and then stored in cryoprotectant at -20°C until staining. For immunohistochemistry, sections were washed three times in 1X PBS for 5 min and permeabilized for 30 min with 0.01% Triton-X in 1X PBS. Sections were next blocked with 10% NGS (Abcam) for 1 h and incubated overnight with primary antibodies at 4°C in 2% NGS. The next day, they were washed three times in 1X PBS before incubating with secondary antibodies for 2 h at room temperature. Sections were washed three times in 1X PBS for 5 min before being mounted on slides with Fluoromount-G (SouthernBiotech).

To validate reduction of Connexin 43 (Cx43) protein in astrocytes transduced with AAVs to express GCaMP-GFP and Cre-RFP, primary antibodies for α-connexin-43 (1:1500, rabbit, Sigma-Aldrich), α-GFP (1:3000, chicken, Abcam), and α-mCherry (1:2000, rat, Thermo Fisher Scientific) in 2% NGS were used. Secondary antibodies include α-rabbit Alexa Fluor 405, α-chicken Alexa Fluor 488, and α-rat Alexa Fluor 555 (all Thermo Fisher Scientific), which were all used at 1:1000 dilution. 60x multi-channel z-stack images were acquired on a CSU-W1 Spinning Disk Confocal (Nikon) using MicroManager from V1 in which AAVs were injected. To quantify loss of Cx43 in RFP⁺ and RFP^- astrocytes, Fiji (ImageJ) was used. Through batch processing, cell maps were created through a semi-automated pipeline to segment astrocytes, with post hoc ROI adjustments for vasculature artifacts. Multi-channel z-stacks were split into 405, 488, and 555 channels, and unstacked into sequential 8-bit z-plane images. For each z-plane, RFP⁺ and RFP^- astrocytes were detected using a Gaussian blur (sigma = 3), thresholding using the Phansalkar method (radius = 1000), and applying ImageJ’s “Analyze Particles” (size > 175 µm², circularity = 0–0.60) to outline ROIs using the wand tool. Corresponding Cx43 images were binarized and the Fiji plugin SynQuant⁵⁸ was used to detect Cx43 puncta number within each RFP⁺ and RFP^- astrocyte in a z-plane’s cell map. Puncta counts were normalized to astrocyte area, and the normalized count from each z-stack was averaged for each slice.

2P image and data analysis

Individual-astrocyte cell maps for time series images were created in Fiji using the following process: For each FOV, an 8-bit z-projection of the time series was created. The z-projection was binarized using the ‘Auto Local Threshold’ feature, using the Niblack method and a radius of 30 or 75, for 16× and 40× images, respectively. Cell maps were drawn on binarized images using a combination of the Lasso and Blow Tool and freehand drawing tool in Fiji, and verified on the z-projected image. Cell maps were also verified against a static indicator of astrocyte morphology when available (EAAT2-tdT⁺ mice for bath-application of LY379268 and Baclofen; GFAP(0.7)-RFP-T2A-iCre in Cx43^floxed mice). To load cell masks into AQuA, regions were saved to the ROI manager and filled in with a color. The regions were projected onto a black image the same size as the original (512 × 512 pixels). The overlay of regions was flattened, converted to an 8-bit image, and saved as a tiff. For the 12.5-min recordings with sequential rounds uncaging glutamate, drift of the slice in X and Y was corrected post-hoc using moco⁵⁹.

AQuA: GCaMP and GluSnFR 2P image sequences were analyzed using AQuA⁸ and custom MATLAB (MATLAB R2018b) and Python (v3.8.18) code. Signal detection thresholds were adjusted for each video to account for differences in noise levels after manually checking for accurate AQuA detection. Cell maps were loaded into AQuA to define cells consistently over multiple time-series featuring the same FOV. For all bath-application experiments, the direction of pia was marked as anterior. For 2P uncaging experiments, the uncaging site was marked as a 3 × 3-pixel landmark.

Bath-application event-based analysis: For Baclofen and t-ACPD Ca²⁺ imaging experiments, Event count per frame was quantified by counting all AQuA-detected events, new or ongoing, in each frame (Fig. 1c). Percent of field active values were calculated by recording the number of active pixels in each frame, as determined by the frame-by-frame footprints of AQuA-detected events. These values were normalized by total number of active pixels in the recording and multiplied by 100. For the Percent of field active dose-response curve (Fig. 1e), the percent of field active values from all frames within the chosen timepoints were averaged by slice. Event propagation was calculated by summing the growing propagation from all cardinal directions, using the AQuA feature propGrowOverall. For dose-response curves for discrete event features (area, duration and propagation [Fig. 1f–h]), all detected Ca²⁺ events within the chosen timepoints were averaged by slice.

The frame the agonist entered the recording chamber was estimated using fluorescence from Alexa Fluor 594 (0.1–2 µM) added to the ACSF reservoir along with agonist. The frame agonist entered the recording chamber was estimated using the maximal curvature method on frames 1–600 of the raw Alexa Fluor 594 fluorescence trace. The maximum curvature method⁶⁰ defines the onset fluorescence changes as the point of maximum curvature during the rising phase of the signal. To identify this point, traces were fit using a modified Boltzmann’s sigmoidal equation:

where a is the difference between the minimum and the maximum fluorescence, b is the inflection point, c is the baseline fluorescence and d is the slope, using a nonlinear least squares algorithm (Levenberg-Marquardt) in MATLAB (Mathworks). Next, the frames of maximum curvature were calculated by setting the fourth derivative of the fitted curve equal to zero and solving for its three solutions. The earliest frame identified out of these three solutions was recorded as the onset frame.

Bath-application ROI-based analysis: Pink Flamindo and GCaMP imaging experiments were analyzed using ROI-based approaches in Fiji. To identify responding cells in Pink Flamindo experiments (Extended Data Fig. 1j), sigmoidal curves were fit to ΔF/F traces using the modified Boltzmann’s sigmoidal equation detailed above. Cells were defined as “responding” if the difference between the minimum and maximum values of the fit curve (a in the Boltzmann’s sigmoidal equation) > baseline noise (3 SD of baseline fluorescence). Responding cells were defined as “increasing” if  and decreasing if  

To identify fluctuations in Pink Flamindo and GCaMP fluorescence (Extended Data Fig. 1k), peaks were detected from ΔF/F traces from individual cells. Peaks were counted if they were 3 SD above the mean baseline fluorescence, had a minimum peak width of 5 frames and a minimum distance of 10 frames between detected peaks. The baseline period was defined as all frames before the frame of agonist entry. For GCaMP, all astrocytes exhibiting ≥ 1 AQuA-detected event during the 10-min recording were run through peak finding. For Pink Flamindo, all detected astrocytes were run through peak finding.

For GCaMP experiments, the frame agonist entered the recording chamber was estimated using the fluorescence from Alexa Fluor 594 (0.1–2 µM) added to the ACSF reservoir along with agonist. Time of agonist entry in the recording chamber was estimated by identifying the first frame Alexa Fluor 594 fluorescence reached ³ 3 SD above baseline mean (frames 1–300); only frames > 375 were considered for evaluation of exceeding the threshold. For Pink Flamindo experiments, dye was not added with agonist to avoid spectral overlap. Time of agonist entry in the recording chamber was estimated by adding 90 frames (the average number of frames for ACSF to travel from the reservoir to the recording chamber) to the frame agonist was added to the reservoir of ACSF

2P uncaging event-based analysis: Individual astrocytes were excluded from analyses (Fig. 2–4, Extended Data Fig. 2–7) if the baseline event rate changed significantly. Changes in baseline event rate for each cell were determined by performing Poisson regression of events in 1-s bins during the period from 90–10s pre-uncaging. Cells with a regression coefficient with p < 0.1 at baseline and with > 5 AQuA-detected events throughout the recording were excluded from all analyses, except for Extended Data Fig. 7d RuBi-glutamate uncaging control. ∆F/F values in raster plots (Fig. 2h and 3c) were calculated using the AQuA output dffMatFilter(:,:,2), the ∆F/F traces from events after removing the contributions from other events in the same location. Cells (Fig. 2h) or local astrocyte networks (Fig. 3c) were sorted based on the onset time of a response following uncaging. A response was defined as the first post-stim peak ≥ threshold (mean baseline DF/F + 3SD), with thresholds calculated by cell or local network using 90–0s before uncaging. For Fig. 3c, the z-score of each local network was calculated on the mean DF/F from AQuA-detected events in the network, using a baseline period of 90–0s before uncaging. For the Sholl-like analysis (Fig. 3h), events were sorted into 50µm bands radiating out from the uncaging site based on the minimum distance between an event and the uncaging site at event onset (using the AQuA output ftsFilter.region.landmarkDist.distPerFrame). 50µm-wide bands began 25mm and end 175µm from uncaging, as events <25mm from the uncaging site likely occur within the stimulated astrocyte and >175mm from the uncaging site can be outside the FOV; see Extended Data Fig. 3i. 90–0s before and 0–150s after uncaging used to calculate change in event number/30s/band. In order to categorize events as propagative or static (Fig. 4d–m and Extended Data Fig. 5b–j, 6 and 7c), the total propagation distance of each event was computed by summing the growing propagation from all cardinal directions, using the AQuA feature propGrowOverall. Events were categorized as propagative if the total propagation distance > 1µm.

Statistics and reproducibility for representative micrographs and spatial heatmaps

Representative micrographs were chosen from experiments repeated with similar results from the following n. Fig. 1b: n = 4 slices, 4 mice; Fig. 2c: n = 72 trials, 12 recordings, 4 slices, 2 mice; Fig. 2f–g: n = 28 astrocytes, 7 slices, 4 mice (note the heterogeneity displayed in Fig. 2h for individual astrocyte responses to NT); Fig. 3b: n = 28 FOV, 7 slices, 4 mice; Fig. 4a: n = 28 FOV, 7 slices, 4 mice; Fig. 4c: n = 15 recordings, 5 mice; Extended Data Fig. 1i: n = 8 slices, 3 mice; Extended Data Fig. 3b: n = 91 FOV, 16 slices, 8 mice; Extended Data Fig. 5a: n = 28 FOV, 7 slices, 4 mice. 

Statistics for Fig. 1–3 and associated Extended Data Figures

All statistical tests used and exact n values can be found for each figure in the corresponding figure legend. Adjustments for multiple comparisons using Bonferroni-Holm correction were implemented using fwer_holmbonf⁶¹. Significance levels defined as the following: ns: p ≥ 0.05, *: p < 0.05, **: p < 0.01, ***: p < 0.001.

Permutation testing: Statistical significance for time-series (t-series) data was computed using permutation testing with custom-written code in MATLAB. 10,000 permutations were run and one- or two- sided p-values for each time point were calculated. P-values were corrected for multiple comparisons using the Benjamini-Yekutieli procedure (implemented using ref. ⁶²) with a False Discovery Rate (FDR) ≤ 0.05.

Data were shuffled/permuted in the following way: To test change in event number/cell (Fig. 1c, Extended Data Fig. 2b and 3g,h), events were shuffled independently for each active cell (≥ 1 AQuA-detected event) in each t-series. For each active cell, events were randomly placed in time bins spanning the duration of the recording (time bins = 60s [Fig. 1c] and 30s [Extended Data Fig. 2b and 3g,h]) and the change in number of events/time bin was calculated as for the experimental data. Permuted changes in event number/cell were averaged across active cells in each t-series and across all t-series to obtain the permuted mean for one round of permutation testing.

To test change in event number/ band (Fig. 3h), permutation tests were run separately for each band and events were shuffled independently for each t-series. For each t-series, events from the tested band were randomly placed in 30 s time bins spanning the duration of the recording, and the change in event number/30 s was calculated as for the experimental data. Permuted changes in event number/30 s were averaged across all t-series to obtain the permuted mean for one round of permutation testing. To test magnitude of change in experimental data versus permuted data, two-sided p-values were calculated as:

(# of times |permuted change| ≥ |experimental change|) +1

# of permutations + 1

For testing increases in ∆F/F (Extended Data Fig. 1d), frames were shuffled independently for each t-series. For each t-series, the average ∆F/F/frame from active regions (≥ 1 AQuA-detected event in either condition [baclofen or t-ACPD]) was calculated, the frame order was shuffled, and the mean ∆F/F/30s was calculated. Permuted mean ∆F/F was averaged across all t-series to obtain the permuted mean for one round of permutation testing. To test magnitude of increases in experimental data versus permuted data, one-sided p-values were calculated as:

(# of times the permuted mean ≥ the experimental mean) +1

# of permutations + 1

Statistics for Fig. 3i­– l, Fig. 4, and associated Extended Data Figures

2P uncaging grid-based ROI analysis: Grid-based regions of interest (ROIs) were determined by dividing the 300 x 300 µm imaging field into a uniform 20 x 20 µm grid (Fig. 3i–l). Each identified Ca²⁺ event was assigned to the ROI in which the centroid of its spatial footprint was located. ROIs with any baseline events were identified as ROIs with ≥1 events in the baseline window 60–0s before uncaging. “Active” ROIs for each NT were identified as ROIs with a ≥ 50% increase in event rate in the window 0–120s after uncaging for that NT, as compared with the rate during the baseline window. Active ROIs were a subset of ROIs with baseline events, as the relative increase in event rate is not defined when there are no baseline events, which results in division by 0. The distance from the uncaging site to each active ROI was determined using the Euclidean distance between the uncaging site, at (0, 0), and the center of each grid ROI (Fig. 3j).

The fraction of overlap (i.e., Jaccard index) O_i between active ROIs for GABA and glutamate were determined for the ith field of view by

where A_GABA,i and A_glutamate,i are the sets of active ROIs for GABA and glutamate, respectively. The overall fraction of overlap O between active ROIs for GABA and glutamate was computed as the mean of the individual O_i  (Fig. 3l).

To determine if the observed fraction of overlap was expected due to chance, a distribution of N = 10,000 surrogate fractions of overlap was computed. The kth surrogate value,  was computed as above, but replacing, for each NT, the set of active ROIs A_NT,i with a new set,   , which was chosen as a random subset of size |A_NT,i| of the set of ROIs with any baseline events for that NT. The p-value for this comparison was estimated⁶³ as

Propagation probability (Fig. 4b): Each Ca²⁺ event was identified as “growing in the depth axis” if the frontier of that event’s spatial footprint extended over time either toward the pia or away from the pia, as determined by the posterior and anterior component of the propGrowOverall metric computed via segmentation by AQuA⁸.

The probability of events growing in the depth axis was computed separately for recordings of GABA and glutamate uncaging within each examined time window. Probabilities were estimated for the baseline window of 60–0s before uncaging, as well as in nonoverlapping 30s bins ranging from 0–150s post-uncaging, by computing the fraction of events that were identified as growing in the depth axis among all events from all recordings within the relevant time window. The change in the probability of events growing in the depth axis was then estimated for each bin as the difference between the fraction of events growing in the depth axis for that bin versus for the baseline period.

To empirically determine the distribution of each of these estimators, we performed this same procedure for estimating the probability of events growing in the depth axis for each NT and time bin on surrogate data generated by hierarchically bootstrapping Ca²⁺ event data, where the hierarchy was sampled cells within sampled recordings (i.e., all events for an individual cell-recording always remained together); this procedure was repeated 10,000 times for each bin. Standard errors were computed as the standard deviation of these empirical distributions.

To determine the probability of observing effects this large under a null hypothesis of no effect of time on the probability of events growing in the depth axis, we computed the distribution of the estimator under an imposed condition in which the overall temporal structure of astrocyte Ca²⁺ events was disrupted. To do this, we performed the same procedure as above for estimating the probability of events growing in the depth axis for each bin, but on surrogate data generated by circularly shifting the timing of each individual cell’s Ca²⁺ events from 90s before to 150s after uncaging by its own independent, uniform random shift between 0s and 240s; this procedure was repeated N = 10,000 times for each bin. As it was unknown whether event propagation would increase or decrease post-uncaging, two-sided p-values were estimated⁶³ as where X denotes the actual observed value of the estimator, and each  is the value of the estimator computed from the kth shifted dataset. These p-values were adjusted across tested time bins and NTs using the Benjamini-Hochberg procedure to obtain q-values, as implemented in statsmodels 0.12.2 (ref ⁶⁴).

Event feature changes (Extended Data Fig. 4a,b): Each Ca²⁺ event is assigned several metrics by AQuA-segmentation⁸, including size (area, perimeter, circMetric [circularity, based on area and perimeter]), amplitude (dffMax), and dynamics (rise19 [rise time], fall91 [fall time], decayTau [decay time constant], width11 [duration]). For each non-propagation metric, the mean metric value among events was computed separately for recordings of GABA and glutamate uncaging for the baseline window 60–0s before uncaging, as well as in nonoverlapping 30s bins from 0–150s post-uncaging. For each bin, the ratio of that bin’s mean metric value to the baseline mean metric value was computed.

AQuA metrics also capture information about events’ directional propagation. Each Ca²⁺ event was identified as “growing” or “shrinking” in each cardinal direction if the frontier of that event’s spatial footprint extended or receded, respectively, over time in that direction, as determined by the components of the propGrowOverall and propShrinkOverall metrics. For each propagation metric, the change in the probability of events growing or shrinking in each axis was computed separately for recordings of GABA and glutamate uncaging within each examined time window, as above in Propagation probability, but using the “growing” or “shrinking” identifiers for each cardinal direction.

To empirically determine the distribution of each of these estimators (i.e., binned post/baseline ratio for non-propagation metrics, binned change in growing/shrinking probability for propagation metrics), we performed the same procedures for computing each metric’s relevant estimators for each NT and time bin outlined above on 10,000 surrogate datasets generated by hierarchically bootstrapping Ca²⁺ event data, as described in Propagation probability. Standard errors were computed as the standard deviation of these empirical distributions.

To determine the probability of observing effects this large under a null hypothesis of no effect of time on the probability of events growing in the depth axis, we computed the distribution of each estimator under 10,000 realizations of an imposed condition in which the overall temporal structure of astrocyte Ca²⁺ events was disrupted by randomly circularly shifting each cell’s Ca²⁺ events, as described in Propagation probability. As it was unknown whether event propagation would increase or decrease post-uncaging, two-sided p-values were estimated using equation (2) above⁶³. These p-values were adjusted across tested time bins and NTs using the Benjamini-Hochberg procedure to obtain q-values, as implemented in statsmodels 0.12.2 (ref. ⁶⁴).

Comparison of in vivo and ex vivo event propagation (Fig. 4d): Events were categorized as propagative or static, as outlined above in the 2P uncaging event-based analysis section. The fraction of propagative events observed in vivo and ex vivo was calculated using baseline events. Ca²⁺ events in in vivo recordings were labeled as “baseline events” if they occurred during periods when the mouse was stationary, as outlined above in the in vivo 2P imaging section. Ca²⁺ events in ex vivo recording were labeled as “baseline events” if they occurred in neighboring astrocytes (i.e. cells not directly stimulated by NT) during the 60–0s before NT uncaging.

To determine the distribution of the two median propagative event fractions empirically, we computed the medians of 10,000 bootstrapped samples of the per-recording fractions for each setting. Standard errors for each statistic were determined from the standard deviations of these empirical distributions.

Computing rate changes for propagative and static events (Fig. 4f, j and Extended Data Fig. 6b–c): The overall rates of propagative and static events for neighboring astrocytes were computed separately for recordings of GABA and glutamate uncaging.

For each event class (i.e., propagative and static events), for each recording, the event rate was computed in each time window as the total number of events from all neighboring cells in that recording in the given time window divided by the duration of that time window. These recording-level rates were computed for the baseline window of 60–0s before uncaging and in nonoverlapping 30s bins ranging from 0–150s post-uncaging. For each recording, the relative rate of propagative and static events was computed for each time bin as the ratio of the event rate for the given event class in that time bin divided by the corresponding event rate in the baseline window. For each time bin, the overall relative rate was estimated as the median of the per-recording relative rates in that time bin.

To determine the distribution of each of these relative rate estimators empirically, we performed this same procedure for estimating relative event rates on surrogate data generated by hierarchically bootstrapping Ca²⁺ event data 10,000 times for each bin (as above in Propagation probability). Standard errors were computed as the standard deviation of these empirical distributions.

To determine the probability of observing effects this large under a null hypothesis of no effect of time post-uncaging on the rate of astrocyte Ca²⁺ events, we computed the distribution of the relative rate estimators under an imposed condition in which the overall temporal structure of astrocyte Ca²⁺ events was disrupted via a random circular shift of the events in each cell, as above in Fig. 4b; this procedure was repeated N = 10,000 times for each bin. Motivated by results in bath application experiments above demonstrating robust aggregate astrocyte Ca²⁺ activity increases in response to agonism of glutamate receptors (Fig. 1h), one-sided p-values were estimated from these permuted datasets, as in equation (1) above. These p-values were adjusted across tested time bins and NTs using the Benjamini-Hochberg procedure to obtain q-values, as implemented in statsmodels 0.12.2 (ref. ⁶⁴).

Determining responding cells based on static and propagative events (Fig. 4h,k and Extended Data Fig. 6e–f): The overall rates of propagative and static events were computed for each neighboring astrocyte, with paired measurements made for recordings of GABA and glutamate uncaging. For each neighboring astrocyte, for each event class (i.e. propagative and static events), the event rate was computed in each time window as the total number of events from that cell in the given time window divided by window’s duration (baseline window: 60–0s before uncaging, response window: 0–120s after NT-uncaging; Extended Data Fig. 5c). Relative event rates were calculated as for Fig. 4f, j and Extended Data Fig. 6b–c above. Cell-recording combinations with zero events of a given type in the baseline window were excluded for computation of relative rates of propagative (GABA: 36 cell-recordings [26.7% of total]; glutamate: 37 [32.2%]) and static (GABA: 0; glutamate: 0) events, as the relative rate would require a division by zero and be undefined in those cases. Astrocytes were identified as “responders” with a particular event type (i.e., static or propagative) to GABA or glutamate if their relative rate of that type of event was ≥ 1.5 for the corresponding NT uncaging recording (Extended Data Fig. 5d). The fraction of astrocytes that were responders was computed for each individual recording, as well as the overall fraction of responders averaged across all recordings for each NT.

To determine the distribution of these overall responder fractions, we performed this same procedure for estimating relative event rates on surrogate data generated by hierarchically bootstrapping Ca²⁺ event data 10,000 times (as above in Propagation probability). Standard errors were computed as the standard deviation of these empirical distributions.

To determine whether there were significant differences between the overall responder fractions for GABA and glutamate, we computed the distribution of the difference between these two fractions under an imposed condition in which there was no systematic difference between GABA and glutamate. To do this, we performed the same procedure as above for estimating the difference between the overall responder fractions for “GABA” and “glutamate”, but on surrogate data generated by, for each cell, swapping the labels for “GABA” and “glutamate” responses from that in the experimental data with probability 1/2; this procedure was repeated 10,000 times. As it was unknown a priori whether GABA or glutamate would have a higher fraction of responder cells, a two-sided p-value was estimated as in equation (2) above.

Decoding NT identity from propagative event responses (Fig. 4i): To quantify the extent to which the observed difference in propagative event responses to uncaged glutamate and GABA enabled reliable identification of NT identity on a trial-by-trial basis, we built a simple classifier that took as input a single value, the relative change in propagative event rate across a FOV in the window 0–120s post-uncaging relative to the window 60–0s pre-uncaging, and classified that FOV as responding to glutamate if the value was ≥ a set threshold, and GABA if the value was < the threshold. To evaluate this classifier’s performance, we built a receiver operating characteristic (ROC) curve by varying the classification threshold across the entire domain of the feature, and at each value of the threshold, computing the empirical true positive rate and false negative rate of the classifier. With the threshold fixed in the ROC analysis, the classifier did not have any remaining free parameters, so did not need to be trained on data and was therefore not a function of any of the data, obviating the need for cross-validation. We computed the area under the ROC curve (AUC) using the trapezoidal rule. To determine the distribution of the observed AUC statistic, we performed this same analysis on 10,000 surrogate datasets generated by bootstrapping (i.e., resampling FOVs with replacement). To determine whether the observed AUC statistic was above 0.5 (indicating completely non-informative decoding) to a degree greater than expected by chance alone, we performed this same analysis on 10,000 surrogate datasets generated by permuting the NT labels.

Determining correlations between GABA and glutamate responses (Fig. 4l): To determine whether individual cells’ responses to GABA and glutamate—as determined in 4h above—were correlated, we computed the Spearman  between the binary paired responses to GABA and glutamate across cells which could be assessed in both conditions (i.e., had > 0 propagating baseline Ca²⁺ events in both recordings) using SciPy 1.6.2 ⁶⁵. To determine the probability of observing a correlation at least this large under a null hypothesis of independence between cells’ responses for GABA and glutamate, we computed Spearman  on surrogate data in which the identities of the cells’ responses to GABA and glutamate were independently permuted; this procedure was repeated 10,000 times. To maintain the ability to identify correlation or anticorrelation, we estimated a two-sided p-value from these surrogate values, as in equation (2).

To complement this analysis, we computed the fraction of overlap (i.e., Jaccard index) between the sets  and  of cells that were responders to GABA and glutamate, respectively:

This statistic is larger when the fraction of overlap between responders for the two neurotransmitters is larger. To determine the probability of observing an overlap at least this large under a null hypothesis of independent responses for GABA and glutamate, we computed this same statistic, but on 10,000 permuted surrogate datasets, as above. To determine significant overlap, we estimated a one-sided p-value from these surrogate values, as in equation (1).

Segregating responding cells based on baseline propagation (Fig. 4m): For each neighboring astrocyte with propagative events during the baseline period of 60–0s pre-uncaging, we computed the fraction of baseline events that were propagative (# propagative baseline events / # all baseline events). Separately for GABA and glutamate, we used the propagative fraction across all given astrocytes to define the threshold fraction of baseline propagative activity, f₅₀, as the 50^th percentile of all observed values; cells with fractions < f₅₀ were said to have “low fraction of propagative events at baseline”, while cells with fractions ≥ f₅₀ were said to have “high fraction of propagative events at baseline” (Extended Data Fig. 5e, top). The fraction of astrocytes that were responders with propagative events to GABA or glutamate were separately estimated from amongst those astrocytes that had low baseline propagation and those that had high baseline propagation, as described above in Determining responding cells based on static and propagative events. Due to the low number of cells in each split group for individual FOVs, the overall average was estimated by pooling all neighboring astrocytes in each group across FOVs.

Similarly for each neighboring astrocyte with baseline propagative events, we computed the rate of all events within the baseline period. Separately for GABA and glutamate, we used the baseline event rate across all neighboring astrocytes to define the threshold baseline event rate, r₅₀, as the 50^th percentile of all observed values; cells with baseline rates < r₅₀ were said to have “low overall baseline event rates”, while cells with fractions ≥ r₅₀ were said to have “high overall baseline event rates” (Extended Data Fig. 5e, bottom). The fraction of astrocytes that were responders with propagative events to GABA or glutamate were separately estimated from amongst those astrocytes that had low overall baseline event rates and those that had high overall baseline event rates, as above.

To determine the distribution of these responder fractions (amongst astrocytes with low and high fraction of propagative events at baseline, or amongst astrocytes with low and high overall baseline event rates), we performed the same procedure for estimating these fractions on surrogate data generated by hierarchically bootstrapping Ca²⁺ event data 10,000 times (as above in Propagation probability). Standard errors were computed as the standard deviation of these empirical distributions.

For each NT, we next sought to determine whether there were significant differences between the fraction of astrocytes that were responders with propagative events amongst cells within the two groupings (i.e., low vs. high fraction of propagative events at baseline; low vs. high overall baseline event rate). Separately for GABA and glutamate, for each group comparison, we computed the difference between the two responder fractions, as well as the distribution of this difference under an imposed condition in which there was no systematic difference in uncaging response between astrocytes in the two groups. To do this, we performed the same procedure as above for estimating responder fractions in the specified groups (e.g., “low fraction of propagative events at baseline” and “high fraction of propagative events at baseline”) as well as the difference between the two, but on surrogate data generated by permuting the group labels; this procedure was repeated 10,000 times. As it was unknown a priori which group in either comparison—low or high baseline propagation, or low or high overall baseline event rate—would have a higher fraction of responder cells, a two-sided p-value was estimated from these surrogate values, as in equation (2).

Simulations to validate characteristics of responder fraction estimates (Extended Data Fig. 5k): Stratifying propagative event responses by the fraction of propagative events in the baseline may induce regression to the mean (RTM) effects, resulting in a bias toward higher observed responsiveness in the low fraction of propagative events at baseline group as compared to the high fraction group. In general, observed effects in differences of repeated measurements stratified by baseline values can arise from a combination of RTM effects and real effects—with the strength of the contribution from RTM depending on the dependency structure and measurement error characteristics in the data—complicating attribution of the observed total effect. To contextualize the observed effect sizes relative to the distribution of effects produced from a pure RTM process, we performed the same procedure as above for estimating responder fractions in the low and high fraction of propagative events at baseline groups, but using surrogate data generated using a random point process model. This model produced simulated event data structured in the same way as the observed dataset: for each cell, the model generated two independent homogeneous Poisson processes, one corresponding to static events and the other corresponding to propagative events. During the simulated baseline period, from 60s to 0s pre-”uncaging”, the rates of these two processes in each cell were set to the observed rate of the corresponding type of event during the veridical baseline period. During the simulated post-”uncaging” period, from 0s to 120s, the rates of these two processes in each cell were determined by multiplying that cell’s baseline rate for the corresponding event type by a response ratio, which was chosen from the empirical distribution of observed post-/pre-uncaging event ratios from amongst all neighboring cells for the given event type. In this way, the simulation modeled the overall response characteristics for propagative events, but in a way that was decoupled from the propagative event fraction in the baseline period.

This simulation procedure was repeated 10,000 times, resulting in a distribution of low-high response fraction differences observed in surrogate data structured in the same way as either the GABA or glutamate uncaging datasets, but with no explicit dependence of cells’ propagative event responses on the baseline propagative event fraction. To summarize the observed effect relative to the effects seen in these simulations, we calculated the fraction of simulations with low-high differences larger than the observed effect.