Unraveling ChR2-driven stochastic Ca2+ dynamics in astrocytes: A call for new interventional paradigms.

Optogenetic targeting of astrocytes provides a robust experimental model to differentially induce Ca2+ signals in astrocytes in vivo. However, a systematic study quantifying the response of optogenetically modified astrocytes to light is yet to be performed. Here, we propose a novel stochastic model of Ca2+ dynamics in astrocytes that incorporates a light sensitive component-channelrhodopsin 2 (ChR2). Utilizing this model, we investigated the effect of different light stimulation paradigms on cells expressing select variants of ChR2 (wild type, ChETA, and ChRET/TC). Results predict that depending on paradigm specification, astrocytes might undergo drastic changes in their basal Ca2+ level and spiking probability. Furthermore, we performed a global sensitivity analysis to assess the effect of variation in parameters pertinent to the shape of the ChR2 photocurrent on astrocytic Ca2+ dynamics. Results suggest that directing variants towards the first open state of the ChR2 photocycle (o1) enhances spiking activity in astrocytes during optical stimulation. Evaluation of the effect of Ca2+ buffering and coupling coefficient in a network of ChR2-expressing astrocytes demonstrated basal level elevations in the stimulated region and propagation of calcium activity to unstimulated cells. Buffering reduced the diffusion range of Ca2+ within the network, thereby limiting propagation and influencing the activity of astrocytes. Collectively, the framework presented in this study provides valuable information for the selection of light stimulation paradigms that elicit desired astrocytic activity using existing ChR2 constructs, as well as aids in the engineering of future application-oriented optogenetic variants.


Introduction
The role of astrocytic calcium signaling in various regulatory mechanisms in the brain is far from being fully understood and is a subject of considerable controversy [see refs [1][2][3][4][5][6]]. These are fueled in part by the limited in vivo reproducibility of in vitro/in situ experimental observations [7][8][9][10], as well as a dearth of cell-specific protocols to induce astrocytic Ca 2+ signaling in vivo in order to delineate their role from ongoing neuronal activity. Sensory and transcranial direct current stimulation techniques [11][12][13] have been used to elicit Ca 2+ changes in astrocytes. However, these methods lack cell specificity due to the concurrent activation of other cell types, including neurons. Studies have also modulated astrocytic Ca 2+ activity in vivo using cell-specific techniques including Ca 2+ uncaging [14,15], chemogenetics [16,17] and optogenetics [18,19]. Ca 2+ uncaging requires invasive site-specific delivery of calcium vehicles, which are rendered inoperative upon depletion via photolysis. Chemogenetics, e.g., designer receptor exclusively activated by designer drugs (DREADDS), offers a platform for controlled, targeted drug delivery; however, they suffer from low temporal resolution. Contrarily, optogenetics is an avant-garde, minimally invasive, and reproducible approach [20,21], providing a platform to genetically target specific cell types with high temporal and spatial precision, which can be employed as a tool to exclusively modulate astrocytic Ca 2+ signaling in vivo.
Despite the recent inception of the field of optogenetics, a wide variety of optogenetic tools have been constructed, among which channelrhodopsin 2 (ChR2) has been one of the most commonly used. Applications of this technique in astrocytes has helped examine the role of these cells in memory enhancement [22], cortical state switching [23], and hyperemic response [18,19]. The biophysical characterization and the response to light stimulation in several ChR2 variants, predominantly in excitable cells, are available in literature [24,25]. More recently, several ChR2 variants have been engineered for enhanced channel conductance (ChETA) [26], increased calcium permeability (CaTCh) [27], and faster recovery kinetics (ChRET/TC) [28]. ChR2 constructs have also been modified to form chimeric variants to regulate responses and facilitate multiwavelength optogenetics in neurons [29]. Contrary to neurons, a holistic approach to quantify the effect of light stimulation on astrocytes has not yet been formulated. Given that light activation of ChR2-enabled astrocytes alters dynamics of intracellular ionic species, mathematical modeling can be of importance in predicting how laser specifications, as well as the biophysical properties of the ChR2 construct, can affect astrocytic calcium signaling. Available theoretical models of Ca 2+ dynamics in astrocytes primarily rely on elevating intracellular IP 3 levels to initiate a Ca 2+ response from the intracellular stores and have not evaluated the effect of direct influx of these ions through transmembrane channels, e.g. ChR2. To achieve predictions with high accuracy, it is also imperative that the model accounts for the stochastic nature of spontaneous calcium oscillations in these cells, which result from the random opening of the inositol trisphosphate receptor (IP 3 R) channels in clusters. Such a model can guide experimentalists in optimizing light stimulation paradigms for existing optogenetic variants to achieve desired Ca 2+ levels in astrocytes, as well as aid in the development of novel application-based constructs targeting these cells.
To this end, we outline a novel stochastic model of astrocytic calcium dynamics with an incorporated optogenetic component-ChR2. Firstly, we quantify and evaluate the effect of different light stimulation paradigms on the Ca 2+ dynamics of single cells expressing three existing ChR2 variants, i.e., wild type (WT), ChETA, and ChRET/TC. For the WT variant, we use the two channel characterizations provided in ref. [30], namely WT 1 and WT 2 . Using the model, we analyze the effect of light stimulation in astrocytes congruent with experimental recording durations (in the order of tens of minutes) [31][32][33], and also gauge the potential effect of long-term optogenetic stimulation (in the order of several hours or beyond) in these cells. Secondly, we quantify the effect of varying stimulation light intensities on a ChR2-incorporated astrocyte with respect to its spiking rate and basal level. Thirdly, to identify key features necessary for the development of prospective ChR2 constructs, we perform a global sensitivity analysis of different parameters of the single cell model to the model output. Lastly, through the incorporation of gap junctions allowing for the diffusion of IP 3 and Ca 2+ , we analyze the effect of local light stimulation on the global Ca 2+ response in a network of astrocytes expressing ChR2.

Materials and methods
The biophysical model outlined in Fig 1 consists of our previously developed model of astrocytic calcium dynamics [34,35], based on the Li-Rinzel simplification of the De Young-Keizer model [36,37], and a 4-state model of ChR2 photocurrent kinetics adapted from refs. [30,38]. We incorporated a 4-state model of ChR2 photocycle, as previous studies have demonstrated its superiority in capturing the dynamics of the ChR2 photocurrent compared to 3-state models [30,39]. The ChR2 model assumes two sets of intra-transitional closed/open states, i.e., dark-adapted (c 1 and o 1 ) and light-adapted (c 2 and o 2 ), in describing the dynamics of a ChR2 photocycle. Light stimulation window (S 0 (t)) is modeled as a pulsed train with unit amplitude and is characterized by a period (T) and a pulse width (δ) (expressed as a percentage of T), indicating the duration within the pulse period for which the light is on.
The model incorporates Ca 2+ influx from the extracellular space through a light-evoked ChR2 flux (J ChR2 ), a capacitive calcium entry (CCE) flux via store-operated Ca 2+ channels (SOC) (J CCE ), and a leak flux representing other transmembrane Ca 2+ channels and exchangers (J in ). It also accounts for the release of Ca 2+ from the endoplasmic reticulum (ER) into the cytosol through IP 3 Rs (J IP 3 R ), a form of Ca 2+ -induced Ca 2+ release (CICR) process, and a leak flux (J Leak ). Replenishment of ER Ca 2+ is done via activation of the sarco(endo)plasmic reticulum Ca 2+ -ATPase (SERCA) pump, while plasma membrane Ca 2+ ATPase (PMCA) pump extrudes cytosolic Ca 2+ into the extracellular space. The dynamics of IP 3 concentration in the soma ([IP 3 ]) is included via Ca 2+ -dependent activation of the Phospholipase C δ1 (PLC δ1 ) pathway. Buffering of Ca 2+ is explicitly accounted for by a fast buffering approximation process [40]. Intercellular signaling in a network of astrocytes is modeled by incorporation of Ca 2+ and IP 3 permeable gap junctions between neighboring cells.

Single cell model
The dynamics of a single ChR2-enabled astrocyte can be summarized by the following system of stochastic differential equations (SDEs): in detail in later sections. G is a diagonal matrix of diffusion rate components of the stochastic processes (noise) associated with state variables. The Brownian motion vector is denoted as dω.
The components of the G matrix are generally determined by the rate parameters of the deterministic part of the SDE, which lead to a state-dependent multiplicative noise in the system [32,44]. In calcium dynamics, the source of stochasticity is mainly attributed to the aggregation of a small number of IP 3 R channels in clusters, which leads to the formation of intracluster Ca 2+ blips (resulting from random opening of a single IP 3 R) or Ca 2+ puffs and sparks (from the concurrent opening of a large number of IP 3 Rs within a cluster). When coupled  Table 1). These processes account for potential sources of stochasticity other than the dynamics of the IP 3 R channels. We further assume that the dynamics of ChR2 gating is only deterministic, i.e., The dynamics of the free cytosolic calcium concentration is given by: where θ = b t K/([Ca 2+ ] c +K) 2 is the buffering factor; b t is the total buffer protein concentration; and K is the buffer rate constant ratio [40]. The efflux of Ca 2+ from the ER into the cytosol via the IP 3 R channel can be described as: The leak of Ca 2+ ions from ER into the cytosol is modeled as: A hill-type kinetic model describing the activity of the SERCA pump is given by: The flux through SOC channels is described using the following equation: Ca 2+ extrusion across the PM via PMCA is given by: PLC δ1 -mediated IP 3 changes in the cell is described as: where X IP 3 denotes the basal rate of IP 3 production in the cell. PLC δ1 activity is described with a Hill's kinetic model as: The dynamics of the fraction of open inactivation IP 3 R gates is given by: where the opening (α h ) and closing (β h ) rates are defined as: where N is the number of IP 3 Rs within a cluster. The total free Ca 2+ in the cell is modeled as: with a zero-mean Gaussian noise component with a constant variance s 2 C o . The open and closed gating dynamics of ChR2 are given by: The existence of ChR2 in open and closed states satisfies the following algebraic condition: The current generated through ChR2 is given by: where G(V m ) is the voltage dependent rectification function determining the shape of ChR2 photocurrent with changing V m . The resultant flux through ChR2 is then derived as Although ChR2 is a non-selective cation channel, in this study, we assume that J ChR2 is solely a calcium flux, as calcium dynamics is one of the most prominent modes of signaling in astrocytes.

Quantification of astrocytic activity
Throughout this manuscript, we calculated the light-evoked mean spiking rate and basal level of [Ca 2+ ] c as a measure of astrocytic activity. For spiking rate calculations, we first removed the trend of simulated calcium traces, i.e., ChR2-induced steady rise in the baseline levels as observed in the traces of Figs 2 and 3. Using a threshold of 0.2 μM, calcium spikes were detected from the baseline-corrected traces and used to compute the mean spiking rate over the entire duration of stimulation. This threshold was chosen to exclude small-amplitude Ca 2+ fluctuations from the calculation of the mean spiking rate. Light-induced elevations in the Ca 2+ basal levels were calculated at the end of the stimulation window. These two measures were compared to study the effect of different light stimulation paradigms, i.e., varying combinations of T and δ values, as well as model parameters, on the response of optogeneticallystimulated astrocytes.

Sensitivity analysis
A global sensitivity analysis was performed to assess the sensitivity of the single cell model output, i.e., Ca 2+ spiking rate and basal level, to variation in parameters determining the magnitude of ChR2 photocurrent and those describing the Ca 2+ buffering kinetics. Each parameter was allowed to vary around its control value within a lower and upper bound (Table 2), and the Latin hypercube sampling (LHS) method with uniform distribution was used to select 1000 random parameter sets to perform the sensitivity analysis [47,48]. The upper and lower bounds were selected such that the range encompasses reported values of parameters for different variants of ChR2. The single cell model was solved for each parameter set, and the partial rank correlation coefficient (PRCC) for each parameter was then computed to determine the magnitude of the parameter influence (positive or negative) on the desired model output. A 95% confidence interval was chosen to determine if the exerted influence is statistically significant.

Astrocytic network model
To study the network-wide response of astrocytes to light stimulation, we incorporated Ca 2+ and IP 3 permeable gap junctions between cells in a network of astrocytes. Fluxes: represent the flow of Ca 2+ and IP 3 , respectively, from astrocyte 'i' in the network to its neighboring cells (indicated by index k). Local and global calcium events are simulated upon focal light stimulation of astrocytes, i.e., a region in the center of the network. Simulations were performed with varying levels of Ca 2+ coupling coefficient, D Ca 2þ , in the presence and absence of Ca 2+ buffering. Calcium responses in the network, i.e., the basal levels and spiking rates, were fitted with a symmetrical 2D Gaussian profile (H(z)) to quantify the peak Ca 2+ -basal/spiking (h max,z ), as well as the magnitude of the spread of Ca 2+ -basal/spiking activity (σ z ) within the network: where z represents calcium basal level or mean spiking rate, h min,z is the minimum value of the z in the network, x 0 and y 0 are the indices of the astrocyte at the center, and x and y are the location of cells away from the origin in both directions. The effect of calcium buffering, varying levels of D Ca 2þ , and the number of cells expressing ChR2 on both the peak and spread of basal/spiking activity was investigated. All simulations were performed in MATLAB using the Local Linearization method described in refs. [46,49] with an integration step size of 0.1 ms. A listing of all parameters and their descriptions is provided in Tables 1 and 2. The codes associated with this study can be downloaded from http://web.eng.fiu.edu/jrieradi/CaAnalysisCode/.

Effect of ChR2 variants on the light-evoked response of astrocytes
Representative traces for the time evolution of state variables of Eq 1 for astrocytes expressing different variants of ChR2 are illustrated for 90 minutes of light stimulation (Figs 2 and S1).

Effect of light stimulation paradigm on single cell model response
Heat maps (Fig 3A) Fig 2). This effect is demonstrated in representative traces of Fig 3B, where depending on the T-δ combination, [Ca 2+ ] c may reach supraphysiological, toxic levels which bring the cell outside of the window for regular spiking (i.e., [Ca 2+ ] c window for IP 3 R-mediated Hopf bifurcation; notice the bottom trace of Fig 3B). Consequently, T-δ combinations that elicit high spiking activity in astrocytes during shortterm stimulation (45 min) will transition into regions with medium and low activity as the stimulus progresses (notice the^symbol in spiking rate heat maps, bottom panel of Fig 3A). Conversely, as seen in the top trace of Fig 3B, combinations with low levels of spiking in shortterm stimulations might transition into regions eliciting high calcium activation, with basal levels within the physiological range, as the stimulation progresses. Thus, to remain within optimal spiking and basal [Ca 2+ ] c levels, one might choose stimulus waveforms based on the desired stimulation duration. It should be noted, however, that the predicted behavior is also a consequence of Ca 2+ buffering capacity of the cell, light intensity, and parameters determining the shape and magnitude of the ChR2 photocurrent (Fig 4).
To quantify the regularity of light-evoked calcium spiking in the astrocyte model, we demonstrate the σ ISI −μ ISI relationship for the inter-spike interval (ISI) values, i.e., the time between subsequent calcium spikes, of simulations in Fig 3A. The average and standard deviation of ISI of the [Ca 2+ ] c traces are represented by μ ISI and σ ISI , respectively. This analysis has been performed for different cell types, including astrocytes, during spontaneous as well as IP 3 -induced calcium oscillations [31,32]. In agreement with these experimental observations, results from our simulations also show a positive and linear correlation between μ ISI and σ ISI values in the astrocytic calcium spiking during light stimulation (Fig 3C). Regions with Ca 2+ levels outside of the predicted range for sustained oscillations, i.e., regions with either low or supraphysiological Ca 2+ basal levels, result in infrequent and irregular spiking with higher μ ISI and σ ISI values.

Sensitivity to model parameters
Simulations predict salient features of light-induced responses in our model astrocyte, i.e., regions with elevated calcium spiking and a steady rise in the Ca 2+ basal level. Results in Fig  4A show the effect of light intensity, or equivalently the amplitude of the waveform, on Ca 2+ response of ChETA-expressing astrocytes for three distinct T-δ combinations. These paradigms correspond to regions with low (△), medium (�) and high (^) spiking activity in the 45-minute stimulation heatmap of Fig 3A, respectively. For all combinations, the Ca 2+ basal level shows an upward trend with increasing light intensity (with varying slopes), indicating a larger Ca 2+ influx through ChR2. Similar trends are observed for spiking rates in the case of low and medium T-δ combinations. For the stimulus paradigm with highest spiking under control conditions (I 0 ), the spiking rate plateaus over a range of intensity values and follows a declining trend upon further increase in the light intensity. This indicates that the associated increase in the basal level leads to a decrease in the firing rate of astrocytes as cytosolic calcium levels exit the region where regular oscillations can occur. These model predictions are similar to experimental observations of monotonic increase in the firing of ChR2-positive neurons stimulated with increasing light intensities [50].
We further performed a global sensitivity analysis to evaluate the effect of parameters determining the dynamics of ChR2 photocycle and buffering capacity on the astrocytic Ca 2+ responses for a low light stimulation paradigm (Fig 4B). Parameters that exerted statistically significant influence on the desired model output are marked with an asterisk. Results indicate differential effects of model parameters on Ca 2+ response of the astrocyte to light stimulation. For instance, increasing the maximum conductance of ChR2 in o 1 state (g 1 ) is positively correlated with both spiking rate and basal levels, as the elevated conductance results in higher magnitudes of Ca 2+ influx through the channel resulting in the elevation of basal levels and further activation of IP 3 Rs. Similarly, e 21 (the rate of transition from o 2 to o 1 ) has a positive correlation with both measures, while the rate of the reverse reaction (e 12 ) is shown to be negatively correlated to the model output. These results suggest that by transitioning the state of ChR2 from o 2 to o 1 , both spiking rate and basal level of Ca 2+ will likely increase. This model prediction can be attributed to the higher conductance of the channel in the o 1 state, the transient region of the ChR2 photocurrent in S2 Fig, compared to o 2 state, the plateau region. Sensitivity analysis of the Ca 2+ buffering parameters demonstrates that with increasing buffering capacity, i.e., increasing total buffer concentration (b t ) or reducing the affinity of buffer proteins to Ca 2+ (K), both basal Ca 2+ level and mean spiking rate of astrocytes will expectedly decrease.

Network-wide astrocytic response to light stimulation
A 10-by-10 network of gap-connected astrocytes was modeled in Fig 5 and S1 Movie. Light stimulation (T = 4s, δ = 45%; D Ca 2þ = 0.1 s -1 ) was applied to a central region (the highlighted box in Fig 5A), and simulations were performed for 45 minutes. Fig 5A heatmaps demonstrate the resulting basal Ca 2+ level in the network after light stimulation, with and without the inclusion of Ca 2+ buffering. Under both conditions, results show increased basal levels in the stimulated region and the propagation of calcium to unstimulated cells. Basal levels reached in focal and distal astrocytes without buffering were drastically higher compared to those with buffering. In the presence of Ca 2+ buffering, only the area in close vicinity of the stimulated region exhibits an increase in the basal level; whereas, in the absence of buffering, even the cells farthest from the stimulation region undergo an increase in calcium. This suggests that buffering reduces the diffusion range of Ca 2+ within the network, thereby limiting propagation. Given the focal and centered stimulation of the network, the distribution of calcium can be suitably quantified using a 2D Gaussian fit (Eq 14). In Fig 5B, the peak basal level, along with the spread from the center (in terms of number of astrocytes), are quantified with varying D Ca 2þ values. In both cases, an intuitive decrease in peak basal level coupled with an increase in the spread is observed as D Ca 2þ increases. Both trends plateaued as D Ca 2þ values reached higher than 0.1 s -1 .
Similar responses are observed in the mean spiking rate of the astrocytic network upon light stimulation (Fig 5C). Inclusion of calcium buffering limited the spiking activity of  Fig 3A for low (△), medium (�), and high (^) activity regions. Data is shown as mean ± std (trials = 5). B) Global sensitivity analysis of the astrocytic Ca 2+ response to variations in ChR2 and Ca 2+ buffering parameters during light stimulation (Δ: T = 2 s, δ = 1% (0.02 s)). Parameters were allowed to vary within a range (Table 2) and 1000 parameter sets were selected using Latin Hypercube Sampling (LHS) with uniform distribution. The partial rank correlation coefficient (PRCC) of each parameter was calculated as a measure of their effect on the basal level (red) and spiking rate (blue) of a ChR2-expressing astrocyte. � denotes statistically significant (p < 0.05) positive or negative influence. astrocytes only in regions immediately surrounding the stimulated region. On the contrary, without calcium buffering, the spiking activity propagated to unstimulated areas and resulted in higher mean spiking rate values. Analysis of network-wide peak spiking rate and spread with buffering in Fig 5D (left panel) indicates a decreasing trend of peak spiking with D Ca 2þ . The propagation of the spiking rate, however, remained almost unchanged over the entire range of D Ca 2þ values examined. Conversely, when buffering was not included (Fig 5D, right  panel), the increase in D Ca 2þ resulted in an increasing trend in the spread of spiking activity while the peak spiking rate remained at high values throughout. Increasing the number of cells expressing ChR2 in the network also resulted in a steady increase in both spiking rate and calcium basal levels in the network (S4 Fig) when buffering was included. Under no buffer condition, the peak basal level showed a steep incline, while the spread of the basal activity remained constant. Both peak spiking rate and spread showed an upward trend with increased expression level, with peak levels reaching a plateau when higher than 50% of the cells expressed ChR2. Collectively, simulations in Fig 5 highlight the role of calcium buffering in limiting the diffusion of free calcium in the network of astrocytes and thereby reducing both basal and spiking levels of calcium in cells. Additionally, comparison of the network results with those of single cells (Fig 3) reveals that the dispersion of calcium into the neighboring cells through gap junctions drastically reduces the supraphysiological values of calcium concentrations predicted in isolated cells.

Discussion
The past few decades have witnessed an upsurge in research on understanding the role of astrocytes in brain function. Aside from providing structural support to neurons [51,52] and regulating the 'excitatory-inhibitory' neurotransmitter balance [53], the exact contribution of these cells, particularly in processes involving their Ca 2+ signaling, is unknown and is highly debated [see refs [1][2][3][4][5][6]]. For instance, astrocytes are thought to be directly involved in neurovascular coupling by sensing neuronal activity and releasing Ca 2+ -mediated vasoactive agents to relax smooth muscle cells of parenchymal arterioles [14,54,55]. However, recent experimental evidence has thrown the validity of this assumption into disarray by questioning astrocytes' involvement, and the extent of their influence, in the onset of the hyperemic response [16,56,57]. Another controversy involving astrocytes is their potential role in modulating neuronal activity via Ca 2+ -mediated release of gliotransmitters. While several groups present evidence for the active involvement of gliotransmission in neuromodulation, others argue that this process is not likely to occur under physiological conditions [4,5,58]. Ca 2+ homeostasis is also dramatically altered in astrocytes under several pathological conditions and during reactive gliosis [59]. These include elevated Ca 2+ levels, as well as intercellular waves in epilepsy [60], Alzheimer's disease [34,61], and spreading depression [62,63]. Despite the involvement of astrocytes in these multi-mechanistic phenomena, differentiating their specific contributions from concurrent neuronal activity has remained one of the most enduring challenges in the field. Thus, selectively targeting astrocytes, and their Ca 2+ signaling, using advanced techniques like optogenetics can provide assistance in resolving the above controversies and help find answers to their exact roles in health and disease. The mathematical modeling framework outlined in our study is a step in this direction in providing a tool for experimentalists to precisely achieve desired astrocytic Ca 2+ levels to examine their role in the abovementioned phenomena.

Salient features of light-evoked calcium signaling in astrocytes
Simulation results predict that calcium dynamics in astrocytes, as seen in experimental studies [33,64,65], can be heavily regulated by light-induced activation of ChR2. Consistent in all simulations performed in this study, upon light activation, astrocytes underwent increases in their basal calcium level and exhibited changes in the spiking activity (Figs 2 and 3). Our results demonstrate that the extent of the rise in calcium (reflecting the higher magnitude of entry through ChR2 compared to the rate of scavenging by PMCA, SERCA, and buffer proteins) is largely dependent on the specifications of the laser stimulus, i.e., T-δ combination and light intensity, as well as parameters determining the shape of the ChR2 photocurrent (Figs 3 and  4). Additionally, whether astrocytes show high or low spiking activity is also contingent upon the duration of light stimulation. As such, T-δ combinations with high spiking activity of single cells in short-term stimulations may transition to regions with medium or low spiking in longer durations, or vice versa (Fig 3B). These results emphasize the importance of choosing an 'ideal' T and δ combination for the desired short-term or long-term astrocytic activity. An inaccurate selection of these combinations could prompt astrocytes to an unphysiological Ca 2+ signaling regime, which might be detrimental for the health of the cells. Whether these model predictions are physiologically accurate needs to be validated against experimental observations for both short and long-term activation of ChR2. When coupled with other astrocytes in a network, however, the dispersion of calcium to neighboring cells dramatically reduced the basal levels reached in the stimulated region (Fig 5A), with values depending on the magnitude of the calcium and IP 3 coupling coefficients in the network. This indicates that experimental design for optogenetic stimulation of a network of astrocytes, e.g., in vivo recordings, cannot be solely based on predictions drawn from single cells, and that network activity depends on ChR2 expression level (S4 Fig).

Engineering of application-based ChR2 variants
Several research groups have engineered ChR2 variants with distinct characteristics, e.g. enhanced conductance, sensitivity, and faster recovery kinetics [26,28] for specific applications in excitable cells. Results of our study can be useful in the development of future ChR2 constructs for eliciting desired activity targeting astrocytes. More specifically, results of our sensitivity analysis (Fig 4B) indicated that the kinetics of ChR2 photocycle significantly affect the Ca 2+ spiking rate and basal levels. Intuitively, directing ChR2 to the open states (o 1 and o 2 ) from the closed states (c 1 and c 2 ) leads to an increase in astrocytic activity in response to light stimulation. For instance, decreasing G d 1 and G d 2 facilitates the existence of ChR2 in the open states as they are negatively correlated to both basal level and spiking rate. Also, an increase in p 2 drives the system to the open state and is positively correlated to both measures. However, less intuitively, for the simulations performed in our study, increased astrocytic activity was achieved when ChR2 resided more in the o 1 state as compared to the o 2 state. This can be observed as an increase in e 21

Model limitations and future directions
ChR2 is a non-selective cation channel, with varying permeability for Na + , K + , Ca 2+ , and H + across different variants [24,25]. Cationic entry has shown to induce membrane depolarization which can activate voltage-gated calcium channels, although their functional role in astrocytes is a subject of ongoing debate [58], and result in further influx of calcium, potentially activating large conductance calcium activated potassium currents [66]. Activation of these channels change the membrane potential of the cell and can thus affect the magnitude of the ChR2 photocurrent, since ChR2 is V m -sensitive. Our model does not include the dynamics of other major ionic species and does not account for the abovementioned dependencies on the membrane potential. We have also not explicitly accounted for other intracellular compartments, e.g. microdomains and mitochondria, involved in the calcium dynamics of astrocytes. In a recent study [67], it was demonstrated that inclusion of microdomains and mitochondria compartments reduced the calcium and IP 3 levels required for the activation of IP 3 Rs in nonexcitable cells and can potentially affect the ChR2-induced calcium signaling in our model. Another limitation of the model is that, in this study, we adopted the Langevin approach outlined in ref. [43] for the implementation of stochasticity in the dynamics of IP 3 R channels and have not accounted for the detailed diffusion of calcium ions within and between IP 3 R clusters. This would require a system of partial differential equations as demonstrated in ref. [32]. In this study, we sought to provide a minimalistic theoretical framework which can readily be employed by researchers for the investigation of light induced Ca 2+ responses in astrocytes. Combination of the presented model with more detailed models as in [42] and [68] where exhaustive geometry and dynamics of various ionic species are accounted for, can enhance our understanding of the intricacies of the behavior of astrocytes and their response to light. while expressing various ChR2 constructs. Different stimulation paradigms ranging from T = 1-5 s and δ = 0-100% of T (trials = 5) were applied from 100 seconds, until the end of the simulation. Each column corresponds to an evaluated ChR2 variant, i.e. WT 1 (left), WT 2 (center) and ChRET/TC (right). Heat maps of Ca 2+ basal level (top panels) and spiking rate (bottom panels) for T-δ combinations are depicted. Scale bar for each heat map was capped to 3 μM and 0.6 spikes/min, respectively. For spike detection, a threshold of 0.2 μM above the basal level was utilized. Results show a similar basal level and spiking rate distribution for WT 1 and WT 2 . However, ChRET/TC shows a smaller region of increased spiking activity coupled with a larger region of T-δ combinations eliciting supraphysiological basal level changes.  8 s)]. Simulations were conducted for 45 minutes while the percentage of the central astrocytes randomly selected to express ChETA was varied (25, 50, 75, 100%) in the presence (left panel) and absence (right panel) of Ca 2+ buffering. A total of 5 trials were conducted for each expression level. A symmetric 2D Gaussian fit was used to quantify the response, i.e. peak and magnitude of the spread from the stimulated region. Top row of plots shows the average and standard deviation of peak basal level and σ basal level as a function of expression. Bottom row of plots shows the average and standard deviation of peak spiking rate and σ spiking rate as a function of expression. For spike classification, a threshold of 0.2 μM above basal level was selected. Increasing the number of cells expressing ChR2 in the stimulation region resulted in a steady increase in both peak basal level and spiking rate in the presence of buffering, coupled with a steady increase in their corresponding σ values. A similar trend can be observed in the absence of buffering; however, after a threshold, the basal level continues to increase, while the spiking rate plateaus.