Purinoreceptors and ectonucleotidases control ATP-induced calcium waveforms and calcium-dependent responses in microglia: Roles of P2 receptors and CD39 in ATP-stimulated microglia

Adenosine triphosphate (ATP) and its metabolites drive microglia migration and cytokine production by activating P2X- and P2Y- class purinergic receptors. Purinergic receptor activation gives rise to diverse intracellular calcium (Ca2+ signals, or waveforms, that differ in amplitude, duration, and frequency. Whether and how these characteristics of diverse waveforms influence microglia function is not well-established. We developed a computational model trained with data from published primary murine microglia studies. We simulate how purinoreceptors influence Ca2+ signaling and migration, as well as, how purinoreceptor expression modifies these processes. Our simulation confirmed that P2 receptors encode the amplitude and duration of the ATP-induced Ca2+ waveforms. Our simulations also implicate CD39, an ectonucleotidase that rapidly degrades ATP, as a regulator of purinergic receptor-induced Ca2+ responses. Namely, it was necessary to account for CD39 metabolism of ATP to align the model’s predicted purinoreceptor responses with published experimental data. In addition, our modeling results indicate that small Ca2+ transients accompany migration, while large and sustained transients are needed for cytokine responses. Lastly, as a proof-of-principal, we predict Ca2+ transients and cell membrane displacements in a BV2 microglia cell line using published P2 receptor mRNA data to illustrate how our computer model may be extrapolated to other microglia subtypes. These findings provide important insights into how differences in purinergic receptor expression influence microglial responses to ATP.

NTPDase of ATP COS-7 [18,24] S.2 Mathematical Models and Parameters Table S2 displays the physical constants used for the model.
These expressions dictate the DAG, ATP, and Ca2+ dependent activities of PKC and PLC. The d terms reflect ATP-dependent and spontaneous (ATP-independent) activation of PLC to capture Ca 2+ transients triggered by P2Y receptor activation versus rest conditions, respectively, as we will describe further below.
These terms are then used to evaluate state models for the GTP-bound Ga protein, activate P2Y receptors, and DAG/IP3 as follows: The first expression reflects that ATP and UTP can be used to trigger the P2Y receptors that are commonly expressed in microglia, e.g. P2Y2 and P2Y6 [12,15] (see main text for further elaboration). The second expression is chosen to be a constant in order reflect ATP-independent activity at resting conditions that culminates in spontaneous Ca 2+ release from s [30]. We also utilized the Ca 2+ -dependent model for PKC activation from [6], though some PKC isoforms may be activated in a Ca 2+ -independent fashion [19] Lastly, the original CC model reflects PKC-dependent inhibition of G αq , whereas more recent studies instead implicate PKC inhibition of PLC-β [19]. Accordingly, we removed the PKC-dependent G αq inhibition step in favor of PKC/PLC-β inhibition (Eqn S1). PLC-β activity may also be influenced by Ca 2+ [19], but we omitted this potential dependence owing to the lack of data to describe the dependency. Corresponding parameters for these equations are provided in Table S3.
-mediated Ca 2+ release from the ER is described by the following equations from : The x terms represent states of the involved in the channel's gating, such as those activating ER Ca 2+ release (see J 1,yk ). The V terms help determine the Ca 2+ -dependent inactivation of the (see J 2,yk ), which results in a negative feedback inhibition motif that gives rise to Ca 2+ oscillations. The net Ca 2+ released by the , J IP3R reflects the competing activation and inactivation terms. Parameters used in these equations are given in Table S4.
This model assumes Ca 2+ release as a deterministic process, though the Ca 2+ transients resemble stochastic spiking in experiment [26]. We note that the parameters provided for the respective CC and YK models used different units (e.g. uM vs nM), therefore conversion factors were used to align the units.
The corresponding sets of equations are now In this model, there are two sets of equations to reproduce the Ca 2+ transients mediated by the IP 3 R activity.
The first set of the Cutherbertson-Chay and Young-Kaizer model equation Table S4: Parameters associated with S3. Values in the parentheses reflect activity at resting (ATP-independent) conditions.

Parameters Values Units
is less ATP sensitive and results in larger spikes as the receptor is activated by IP3. The second set of the model equations generates ATP-independent Ca 2+ oscillations [30], the amplitude of which is smaller than the oscillations generated by the ATP-dependent model (the first set).
Although there is no direct evidence that IP 3 maintains the Ca 2+ baseline transients, it is fair to say that the IP 3 -mediated pathway is not amplified rather than terminated or silenced due to the absence of ATP. Therefore, instead of developing a new set of mathematical expressions to mimic spontaneous Ca 2+ baseline, we again adapted the CC and YK integrated model with slightly tuned parameters to mimic the randomness of the baseline transients but the mathematical expression remains a deterministic model.
The input and fitted parameters are reflected in Tables S3 and S4) to best reproduce experimentally-measured intracellular Ca 2+ transients in microglia.
The difference between the ATP-dependent and independent equations are marked by asterisk (*) in Eqn S3. The [AT P ] cc in the ATP-dependent mechanism is the trigger for the subsequent pathways, whose magnitude is based on both the concentration of ATP and UTP, since the majority of P 2Y receptor consists of P2Y2 and P2Y6 [12,15]. On the other hand, the ATP-independent mechanism is simply expressed as a constant variable that triggers consistent oscillatory waveform of Ca 2+ in the system to mimic the baseline in microglia [30]. Unless it is listed, the ATP-independent Ca 2+ baseline transients were generated by the identical model with the set of adjusted parameters listed in Table S3 and S4.
The YK model also includes feedback inhibition of IP3R by Ca 2+ [31]. Although PLC-β is calcium-dependent, it is omitted in the model for simplicity [19].
To integrate two models, the calculation performs the unit conversion for [IP 3 ], [Ca 2+ ] i and [Ca 2+ ] ER as shown Eqn S5.

S.2.2 Markov State Modeling based P2X receptor kinetic model
The following parameters were refit from [4] (P2X-only) to reflect simultaneous contributions from P2X and P2Y-class receptors. P 2X4 The mathematical expression is available in the previous work [4].

S.2.3 Homeostasis Equations in Microglia
Ca 2+ -handling in the ER lumen ([Ca 2+ ] ER ) Our model for Ca 2+ handling in the ER(Eqn S7) reflects of the Ca 2+ release via SERCA (J ERtoCyt,SERCA ), buffer by calsequestrin or other ER-resident Ca 2+ binding proteins (R Ca 2+ ·S ), and Ca 2+ release via IP 3 receptors (J IP3R,AT P −dependent and J IP3R,AT P −independent ). The term denoted 'ATP-independent' represents baseline activity for resting microglia. The leak term (J ERtoCy,Leak ) represents residual Ca 2+ from the ER to offset SERCA Ca 2+ uptake. Detailed expressions for these contributions are found in our previous work [4].

Homeostasis in Cytoplasm domain ([Ca 2+
] i ) Our model for Ca 2+ handling in the cytoplasm domain (Eqn S8) consists of Ca 2+ uptake via SERCA (J CytoER,SERCA ), buffering by the calmodulin-calcineurin complex (R Ca·F ,R Ca·B , and R Ca·CaM −CN ), and Ca 2+ release via IP 3 receptors (J IP3R,AT P −dependent and J IP3R,AT P −independent ). Plasma membrane contributions to the cytosolic Ca 2+ include the sodium/calcium exchanger (NCX) activity (J N CX ) and P2X receptors (J P 2X7 +J P 2X4 ) introduced in [4]. The leak terms (J ERtoCy,Leak and J ExtoCy,Leak ) ensure the system remains in steady state at rest. Detailed expressions for these contributions are found in our previous work [4].
Parameters for our model are provided in Tables S7-S10.    Table S10: Parameters for CaM-dependent activation of calcineurin. The corresponding equations were adapted from [2,11] and listed in the previous work [4].  Table S11: Parameters Ca 2+ buffer interaction calculations. The equations associated with the rest of the buffers were adapted from [25] and listed in the previous work [4].

Parameters Values Units
Bmax,

S.2.4 Miscellaneous signal transduction
NFAT The NFAT model used in our study was adapted from Cooling et al [5] and implemented to link NFAT activation to TNFα synthesis [33]. The equations listed in the previous work [4] were utilized with the parameters (Table S12) fitted in this study. Table S12: Parameters for NFAT cycle calculations. The corresponding equations were adapted from [5] and listed in the previous work [4].

S.2.5 TNFα
According to a report from Barbera et al [1], matured TNFα is present inside the cell and can be exocytosed following the activation of ionotropic P2 receptors. Our model for this process is provide in Eqn S9. Table S13: Parameters for phosphorylation of p38 calculations. The equations associated with p38 were adapted from [29] and listed in the previous work [4].

Parameters Values Units
[pp38] total 100 k b,pp38 8.51 Ohsawa et al [22,23] determined that P 2Y 12 activation drives migration following ADP treatment (a product of ATP degradation by CD39). [15,21]. Their latter work determined the augmentation of migration via ionotropic P2 receptor stimulation that promotes Ca 2+ entry into the cell. Our model for chemotaxis includes both P 2Y 12 activation and the elevation of intracellular Ca 2+ from P 2X4 and P 2X7 receptors (Eqn S10).
As shown by Ohsawa et al [22,23], PI3K activation and Akt phosphorylation are proportional to P 2Y 12 activity. This process proceeds through the activation of Gi/o. (Eqn S11).
(S11) Chemotaxis in our model scales with the active states of pAkt [22] and CaM [8] :

S.2.7 Degradation of ATP by NTPDase1
Ectonucleoside triphosphate diphosphohydrolase-1 (E-NTPDase), also known as CD39, is a plasma-membrane protein that plays an important role in microglial migration by balancing ATP and adenosine molecules [9]. We have adapted the kinetic model of ATP degradation into ADP and AMP by CD39 introduced in the [18] to determine the ATP and ADP available for the purinergic receptors in our microglial models: (S13) The parameters for these equations controlling nucleotide availability are given in Table S16.  Figure S1: Predicted (black) versus experimentally measured (blue) Ca 2+ responses following 100 µM ATP. [14] .

Figure S2:
A) Schematic for Ca 2+ waveforms generated by P 2X4 and P 2Y 2 in response to 1mM ATP for 10 minutes, which does not include the degradation of ATP by CD39. B) Comparison of predicted (black: moving average with a window size of 8 reported along the left y-axis.) and experimentally-measured [14] (blue-dashed) Ca 2+ transients reported along the right y-axis . Figure S3: A) Intracellular Ca 2+ transients and their corresponding B) oscillation of active G αq with respect to k g,cc that controls the activation of G αq . Figure S4: Comparison between Ca 2+ transients induced by ER Ca 2+ release via IP 3 -mediated pathway at 100 µM and 1 mM ATP concentrations in cytoplasm (A) and ER lumen (B). The WT microglia model was used for this prediction. The faded lines denote the contribution by P 2Y receptor activation that results in ER Ca 2+ release to the cytosolic domain. The data demonstrate the relationship between cytosolic and ER Ca 2+ transients, which suggest that roughly 43.7% of Ca 2+ is drawn from the ER at low ATP vs. 33.3% at high ATP.  Predictions of dephosphorylated NFAT in nucleus over time (A) and TNFα mRNA with various computation configurations and comparison to the previously developed model [4] with respect to amplitude of stimulation (5). All simulations in B) were performed for 5 minutes. The plot is in the unit of scale, whose basis is the maximum increment measured by the current model. Figure S9: Simulation of Ca 2+ transients in BV2 with 100 µM and 1 mM ATP and UTP. Primary microglia with low (A) and high (B) ATP/UTP concentration whereas BV2 cells with low (C) and high (D) concentration of the stimulant. In this case, UTP selectively actives only P2Y2 receptors. [27].