Single cell model

The dynamics of a single ChR2-enabled astrocyte can be summarized by the following system of stochastic differential equations (SDEs):

Table 1

Model parameters.

	Value		Unit	Description	Source
IP3 Dynamics
vδ	0.15		μM/s	Maximum rate of IP3 production (PLCδ1)	[35]
KδCa	0.55		μM	Half saturation constant of Ca2+ (PLCδ1)	[35]
$$	1.25		s-1	IP3 degradation rate	[35]
$$	0.14		μM/s	Basal level of cytosolic IP3 production	[35]
Ca2+ Dynamics
xCCE	0.01		μM/s	Maximum CCE influx	[35]
hCCE	10		μM	Half-inactivation constant for CCE influx	[35]
α1	0.19		~	Volume ratio between ER and cytosol	[35]
vSERCA	0.90		μM/s	Maximum rate constant of SERCA pump	[35]
Kp	0.10		μM	Half-maximal activation of the SERCA pump	[35]
d1	0.13		μM	Dissociation constant for IP3 (IP3R)	[35]
d5	0.08		μM	Ca2+ activation constant (IP3R)	[35]
v1	6		s-1	Ligand-operated IP3R channel flux constant	[35]
v2	0.11		s-1	Ca2+ passive leakage flux constant	[35]
kout	0.50		s-1	Rate constant of Ca2+ extrusion	[35]
λ	1		~	Time scaling factor	[35]
ε	0.01		~	Ratio of PM to ER membrane surface area	[35]
Jin	0.04		μM/s	Passive leakage	[35]
Vm	-70		mV	Membrane voltage	Model assumption
G(Vm)	$$		~	Voltage dependent rectification function	[38]
volcyt	10-12		L	Volume of the cytosol, assuming spherical cell	[41]
Am	4.83×10-6		cm2	Surface area of the astrocyte membrane (calculated using volcyt and assuming spherical cell shape)	Model assumption
F	9.65×104		C/mol	Faraday’s constant
$$	2		~	Valence of Ca2+ ions
bt	200		μM	Total buffer protein concentration	[42]
K	20		μM	Buffer rate constant ratio	[42]
N	20		~	Number of IP3Rs within a cluster	[43]
Gating Parameters
a	0.20		(μMs)-1	Rate constant for Ca2+ binding in IP3 inhibitory site	[35]
d2	1.05		μM	Dissociation constant for Ca2+ inhibition (IP3R)	[35]
d3	0.94		μM	Dissociation constant for IP3 (IP3R)	[35]
Network Dynamics
$$	1		s-1	Coupling coefficient of IP3	Model estimate
$$	0.1		s-1	Coupling coefficient of Ca2+	Model estimate
Wiener Processes
$$	0.02		s-1/2	Variance of Wiener process of [IP3]	[35]
$$	0.01		s-1/2	Variance of Wiener process of [Ca2+]c	[35]
$$	0.01		s-1/2	Variance of Wiener process of [Co]	[35]
ChR2 Parameters
p1	ChRWT1	0.06	ms-1	Maximum excitation rate of c1	[26,28,30]
	ChRWT2	0.12
	ChETA	0.07
	ChRET/TC	0.13
$$	ChRWT1	0.46	ms-1	Rate constant for the o1 to c1 transition	[26,28,30]
	ChRWT2	0.01
	ChETA	0.01
	ChRET/TC	0.01
e12	ChRWT1	0.20	ms-1	Rate constant for the o1 to o2 transition	[26,28,30]
	ChRWT2	4.38
	ChETA	10.51
	ChRET/TC	16.11
e21	ChRWT1	0.01	ms-1	Rate constant for the o2 to o1 transition	[26,28,30]
	ChRWT2	1.60
	ChETA	0.01
	ChRET/TC	1.09
p2	ChRWT1	0.06	ms-1	Maximum excitation rate of c2	[26,28,30]
	ChRWT2	0.01
	ChETA	0.06
	ChRET/TC	0.02
$$	ChRWT1	0.07	ms-1	Rate constant for the o2 to c2 transition	[26,28,30]
	ChRWT2	0.12
	ChETA	0.15
	ChRET/TC	0.13
Gr	ChRWT1	9.35×10-5	ms-1	Recovery rate of the c1 state after light pulse is turned off	[26,28,30]
	ChRWT2	9.35×10-5
	ChETA	1×10-3
	ChRET/TC	3.85×10-4
τChR2	ChRWT1	6.32	ms-1	Activation time of the ChR2 ion channel	[26,28,30]
	ChRWT2	0.50
	ChETA	1.59
	ChRET/TC	0.36
g1	ChRWT1	0.11	mS/cm2	Maximum conductance of the ChR2 ion channel in the o1 state	[26,28,30]
	ChRWT2	0.10
	ChETA	0.88
	ChRET/TC	0.56
γ	ChRWT1	0.03	~	Ratio of maximum conductance of the ChR2 ion channel in the o2 and o1 state $$	[26,28,30]
	ChRWT2	0.02
	ChETA	0.01
	ChRET/TC	0.02
EChR2	All variants	0	mV	Reversal potential of ChR2	[38]

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₃R channels in clusters, which leads to the formation of intra-cluster Ca²⁺ blips (resulting from random opening of a single IP₃R) or Ca²⁺ puffs and sparks (from the concurrent opening of a large number of IP₃Rs within a cluster). When coupled with other clusters, synchronized Ca²⁺ puffs can originate global Ca²⁺ events in the cell [32]. Several schemes have been developed to model the stochastic dynamics of IP₃Rs in different cell types, ranging from models incorporating detailed Markov Chain processes, to those adapting Fokker-Planck approximations (or their equivalent Langevin equations) of the Markov approaches [refer to ref. [45] for a review of the deterministic and stochastic models of IP₃Rs].

In this study, we adopted the Langevin approach outlined in ref. [43] for the dynamics of open IP₃R inactivation gates in Eq 1, where the associated noise is a Gaussian white noise with zero mean and the variance ($$ {\mathrm{\sigma }}_{\mathrm{h}}^2$$) depending on the levels of [Ca²⁺]_c, [IP₃], h, and the number of IP₃Rs in a cluster. We have previously estimated the noise level in the equations for [IP₃], [Ca²⁺]_c, and [C_o] in astrocytes as additive, uncorrelated Wiener processes with zero mean and a constant variance ($$ {\mathrm{\sigma }}_{\mathrm{i}}^2;\mathrm{i}=[{\mathrm{Ca}}_{\mathrm{c}}^{2+},{\mathrm{IP}}_3,{\mathrm{C}}_{\mathrm{o}}]$$) [35] using the local linearization (LL) filter [46] (Table 1). These processes account for potential sources of stochasticity other than the dynamics of the IP₃R channels. We further assume that the dynamics of ChR2 gating is only deterministic, i.e., $$ {\mathrm{\sigma }}_{{\mathrm{o}}_1}={\mathrm{\sigma }}_{{\mathrm{o}}_2}={\mathrm{\sigma }}_{{\mathrm{c}}_2}={\mathrm{\sigma }}_{\mathrm{s}}=0$$.

The dynamics of the free cytosolic calcium concentration is given by:

The efflux of Ca²⁺ from the ER into the cytosol via the IP₃R channel can be described as:

The leak of Ca²⁺ 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²⁺ extrusion across the PM via PMCA is given by:

PLC_δ1-mediated IP₃ changes in the cell is described as:

The dynamics of the fraction of open inactivation IP₃R gates is given by:

The total free Ca²⁺ in the cell is modeled as:

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:

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²⁺]_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²⁺ fluctuations from the calculation of the mean spiking rate. Light-induced elevations in the Ca²⁺ 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 optogenetically-stimulated astrocytes.

Fig 2

Response of ChR2 variants to light stimulation.

Representative traces (90 minutes in duration) of IP₃ concentration ([IP₃]), cytosolic calcium concentration ([Ca²⁺]_c), fraction of open inactivation IP₃R gates (h), and total calcium concentration ([C_o]) for an astrocyte expressing various ChR2 variants (wild type 1 (WT₁), wild type 2 (WT₂), ChETA and ChRET/TC) in response to light stimulation (T = 2 s, δ = 20% (0.4 s), unit amplitude). The blue horizontal solid line indicates the stimulation window.

Fig 3

Response of a ChETA-expressing astrocyte to various light stimulation paradigms.

Simulations were conducted to evaluate astrocytic Ca²⁺ response to different paradigms ranging from T = 1–5 s and δ = 0–100% of T (trials = 5). A) Heat maps of Ca²⁺ basal level (top panel) and spiking rate (bottom panel) for the T-δ combinations assessed. Each column depicts basal level and spiking rate heat maps up to a certain point, i.e. 45, 90, and 270 minutes, to determine the cell response as light stimulation progresses from short-term to long-term. Scale bars for heat maps of basal level and spiking rate were capped to 3 μM and 0.6 spikes/min, respectively. The △ and ◊ symbols correspond to T-δ combinations that result in the traces of panel B. B) Representative [Ca²⁺]_c traces corresponding to T-δ combinations highlighted in A [Top trace (△): T = 2 s, δ = 1% (0.02 s); Bottom trace (◊): T = 4 s, δ = 45% (1.8 s)]. The inset of each trace highlights a 15-minute section to show the detected Ca²⁺ spikes. For spike detection, a threshold of 0.2 μM above the basal level was utilized. Each vertical solid grey line denotes the time points corresponding to heatmaps of panel A. C) Average (μ_ISI) vs. standard deviation of ISI (σ_ISI) where each point corresponds to a single simulated 270-min [Ca²⁺]_c trace. Solid line shows the linear fit between μ_ISI and σ_ISI values. Throughout the manuscript, △, ○, and ◊ represent paradigms resulting in low, medium, and high Ca²⁺ spiking during short-term (45 min) stimulation of astrocytes, respectively.

Sensitivity analysis

A global sensitivity analysis was performed to assess the sensitivity of the single cell model output, i.e., Ca²⁺ spiking rate and basal level, to variation in parameters determining the magnitude of ChR2 photocurrent and those describing the Ca²⁺ 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.

Table 2

Range of parameters for global sensitivity analysis.

Parameter	Range	Unit
ChR2 parameters
p1	51.28−1.50×102	ms-1
$$	8.16−5.47×102	ms-1
e12	163.52−1.93×104	ms-1
e21	4.00−1.92×103	ms-1
p2	10.00−7.69×101	ms-1
$$	56.32−1.81×102	ms-1
Gr	7.48×10-2−1.20	ms-1
τChR2	2.89×10-4−7.6×10-3	ms-1
γ	1.13×10-2−3.66×10-2	~
g1	7.84×10-2−1.05	mS/cm2
Buffering Parameters
bt	160−240	μM
K	16−24	μM

Astrocytic network model

To study the network-wide response of astrocytes to light stimulation, we incorporated Ca²⁺ and IP₃ permeable gap junctions between cells in a network of astrocytes. Fluxes:

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). Simulations reveal a progressive increase in the spiking activity of cytosolic calcium in all variants with the duration of stimulus. Apparent in the traces is a lag between the onset of light stimulation and regions with more regular spiking activity, likely reflecting the time required for [Ca²⁺]_c to reach sufficient levels to activate IP₃R channels. Traces show a steady increase in the total calcium content of the cell, and conversely, a decreasing trend in the fraction of open IP₃R inactivation gates (h). The [IP₃] levels, however, show a rather steady profile for WT and ChETA variants, and predict only a slight increase towards the end of the stimulation window for ChRET/TC. This increase can be attributed to the elevated basal [Ca²⁺]_c in ChRET/TC expressing astrocytes, activating the production of IP₃ via the PLC_δ1 pathway ([Ca²⁺] at half- maximal activation of PLC_δ1 is ~0.5 μM, Table 1). Comparison of results across variants suggests that light-induced changes are steepest in astrocytes expressing ChRET/TC. These model predictions are in line with the higher magnitude of the ChR2 flux in ChRET/TC compared to other variants in both transient and plateau regions of the photocurrent during the light stimulus (S2 Fig). Furthermore, results suggest that for the simulation time in this figure, the mechanism of light-evoked increased spiking in cytosolic calcium, predominantly in WT or ChETA enabled astrocytes, is mainly through the activation of IP₃Rs by increases in [Ca²⁺]_c rather than elevations in [IP₃].

Effect of light stimulation paradigm on single cell model response

Heat maps (Fig 3A) show the effect of various light stimulation paradigms on the spiking rate and Ca²⁺ basal level in ChETA-expressing astrocytes for different stimulation times, i.e., 45, 90, and 270 minutes (for other ChR2 variants refer to S3 Fig). Simulations show that as the stimulus duration increases, the Ca²⁺ basal level rises steadily due to a net Ca²⁺ entry into the cytosol (note the upward trend of [C_o] in the traces of Fig 2). This effect is demonstrated in representative traces of Fig 3B, where depending on the T-δ combination, [Ca²⁺]_c may reach supraphysiological, toxic levels which bring the cell outside of the window for regular spiking (i.e., [Ca²⁺]_c window for IP₃R-mediated Hopf bifurcation; notice the bottom trace of Fig 3B). Consequently, T-δ combinations that elicit high spiking activity in astrocytes during short-term 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 short-term 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²⁺]_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²⁺ buffering capacity of the cell, light intensity, and parameters determining the shape and magnitude of the ChR2 photocurrent (Fig 4).

Fig 4

Model sensitivity to light intensity, Ca²⁺ buffering, and ChR2 parameters.

A) Effect of different light intensities and stimulation paradigms on the basal level and mean spiking rate of a ChETA-expressing astrocyte. Light intensity was varied as a fraction of the control value (I₀). Stimulation paradigms were selected based on the 45-min spiking rate heat map in 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²⁺ response to variations in ChR2 and Ca²⁺ 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.

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²⁺]_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₃-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²⁺ levels outside of the predicted range for sustained oscillations, i.e., regions with either low or supraphysiological Ca²⁺ 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²⁺ basal level. Results in Fig 4A show the effect of light intensity, or equivalently the amplitude of the waveform, on Ca²⁺ 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²⁺ basal level shows an upward trend with increasing light intensity (with varying slopes), indicating a larger Ca²⁺ 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₀), 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²⁺ 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²⁺ response of the astrocyte to light stimulation. For instance, increasing the maximum conductance of ChR2 in o₁ state (g₁) is positively correlated with both spiking rate and basal levels, as the elevated conductance results in higher magnitudes of Ca²⁺ influx through the channel resulting in the elevation of basal levels and further activation of IP₃Rs. Similarly, e₂₁ (the rate of transition from o₂ to o₁) has a positive correlation with both measures, while the rate of the reverse reaction (e₁₂) is shown to be negatively correlated to the model output. These results suggest that by transitioning the state of ChR2 from o₂ to o₁, both spiking rate and basal level of Ca²⁺ will likely increase. This model prediction can be attributed to the higher conductance of the channel in the o₁ state, the transient region of the ChR2 photocurrent in S2 Fig, compared to o₂ state, the plateau region. Sensitivity analysis of the Ca²⁺ 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²⁺ (K), both basal Ca²⁺ 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%; $$ {\mathrm{D}}_{{\mathrm{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²⁺ level in the network after light stimulation, with and without the inclusion of Ca²⁺ 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²⁺ 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²⁺ 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 $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ values. In both cases, an intuitive decrease in peak basal level coupled with an increase in the spread is observed as $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ increases. Both trends plateaued as $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ values reached higher than 0.1 s^-1.

Fig 5

Network-wide astrocytic Ca²⁺ response to light stimulation.

A 10-by-10 network of astrocytes was employed to analyze the response of cells when the central 4-by-4 astrocytes (white square in heatmaps of panels A and C) were stimulated (◊: T = 4 s, δ = 45% (1.8 s), 45 min). See the network organization schematic in Fig 1. Simulations were conducted in the presence (left panel) and absence (right panel) of Ca²⁺ buffering with varying gap junctional Ca²⁺ coupling coefficient ($$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$). A symmetric 2D Gaussian fit was utilized to quantify the response, i.e., peak and magnitude of the spread from the stimulated region. A) Heat maps of network-wide light stimulation-induced Ca²⁺ basal levels. B) Plots of the peak basal level and σ_basal-level obtained from the Gaussian fit with varying $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ values. Vertical dashed grey line denotes the $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ value used to generate the heat maps in panels A and C. C) Heat maps of network-wide spiking rate response corresponding to basal levels in panel A. For spike classification, a threshold of 0.2 μM above the basal level was selected. D) Plots for the peak spiking rate and σ_{spiking rate} with varying $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ values.

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 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 $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$. The propagation of the spiking rate, however, remained almost unchanged over the entire range of $$ {\mathrm{D}}_{{\mathrm{Ca}}^{2+}}$$ values examined. Conversely, when buffering was not included (Fig 5D, right panel), the increase in $$ {\mathrm{D}}_{{\mathrm{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.

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²⁺ 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₃ 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²⁺ spiking rate and basal levels. Intuitively, directing ChR2 to the open states (o₁ and o₂) from the closed states (c₁ and c₂) leads to an increase in astrocytic activity in response to light stimulation. For instance, decreasing $$ {\mathrm{G}}_{{\mathrm{d}}_1}$$ and $$ {\mathrm{G}}_{{\mathrm{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₂ 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₁ state as compared to the o₂ state. This can be observed as an increase in e₂₁ and a decrease in e₁₂ leading to the existence of ChR2 in o₁ (see Fig 1 ChR2 photocycle). The shape of the photocurrents in S2 Fig also confirms that the channel has the highest flux in the o₁ state (the transient phase). Collectively, these results suggest that tailoring new ChR2 constructs such that the photocurrent is directed mainly towards the o₁ state can enhance the astrocytic activity. The same analysis can be performed under different stimulation paradigms and for varying durations of stimulus.

Model limitations and future directions

ChR2 is a non-selective cation channel, with varying permeability for Na⁺, K⁺, Ca²⁺, 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₃ levels required for the activation of IP₃Rs in non-excitable 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₃R channels and have not accounted for the detailed diffusion of calcium ions within and between IP₃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²⁺ 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.