Mathematical modelling of human P2X-mediated plasma membrane electrophysiology and calcium dynamics in microglia

Regulation of cytosolic calcium (Ca2+) dynamics is fundamental to microglial function. Temporal and spatial Ca2+ fluxes are induced from a complicated signal transduction pathway linked to brain ionic homeostasis. In this paper, we develop a novel biophysical model of Ca2+ and sodium (Na+) dynamics in human microglia and evaluate the contribution of purinergic receptors (P2XRs) to both intracellular Ca2+ and Na+ levels in response to agonist/ATP binding. This is the first comprehensive model that integrates P2XRs to predict intricate Ca2+ and Na+ transient responses in microglia. Specifically, a novel compact biophysical model is proposed for the capture of whole-cell patch-clamp currents associated with P2X4 and P2X7 receptors, which is composed of only four state variables. The entire model shows that intricate intracellular ion dynamics arise from the coupled interaction between P2X4 and P2X7 receptors, the Na+/Ca2+ exchanger (NCX), Ca2+ extrusion by the plasma membrane Ca2+ ATPase (PMCA), and Ca2+ and Na+ leak channels. Both P2XRs are modelled as two separate adenosine triphosphate (ATP) gated Ca2+ and Na+ conductance channels, where the stoichiometry is the removal of one Ca2+ for the hydrolysis of one ATP molecule. Two unique sets of model parameters were determined using an evolutionary algorithm to optimise fitting to experimental data for each of the receptors. This allows the proposed model to capture both human P2X7 and P2X4 data (hP2X7 and hP2X4). The model architecture enables a high degree of simplicity, accuracy and predictability of Ca2+ and Na+ dynamics thus providing quantitative insights into different behaviours of intracellular Na+ and Ca2+ which will guide future experimental research. Understanding the interactions between these receptors and other membrane-bound transporters provides a step forward in resolving the qualitative link between purinergic receptors and microglial physiology and their contribution to brain pathology.


Introduction
Recently, microglia have attracted wide attention owing to their ability to undergo a variety of morphological configurations in health and disease [1][2][3]. Motile microglial cells play a key defensive role in the central nervous system (CNS) through their role in clearing pathogens and active maintenance of neurons and synapses. Microglia appearance can range from ramified to amoeboid-like cells dependent on their microenvironment [4,5], however this two state model does not capture the full range of microglia heterogenity [6]. Restructuring of the actin cytoskeleton for directed motility in microglia consists of a complex molecular cascade that involves several membrane-coupled receptors [5]. These receptors enable the microglia to detect subtle concentration changes in their environment in order to begin either process extension or whole-cell chemotaxis. Three major families of purinergic receptors are found in human, rat and mouse microglia. These include adenosine receptors (A 1 , A 2A , A 2B and A 3 ), P2X receptors (P2XRs) (P2X 1 , P2X 4 and P2X 7 ) and P2Y receptors (P2Y 2 , P2Y 6 and P2Y 12-14 ) [7]. Evidence to date suggests that P2X 4, P2X 7 and P2Y 12 are the most important receptors for microglial directed motility both in vitro and in vivo, where spontaneous calcium (Ca 2+ ) transients play a key role [8][9][10][11][12]. P2XRs are plasma membrane trimeric assemblies, which bind ATP for ionic permeation through the plasma membrane. One of the important elements in controlling intracellular calcium ([Ca i 2+ ]) is ATP and its activation of purinergic P2XRs in eukaryotic cells [13,14]. Upon adenosine triphosphate (ATP) binding, a channel pore opens for Ca 2+ and other ionic currents to translocate into the cell. Electrophysiology recordings are used to study the P2XR biophysical properties in response to ATP stimualtion. The main goal of this paper is to develop a biologically faithful model of ATP-triggered P2XR Ca 2+ and sodium (Na + ) signalling transduction and plasma membrane electrophysiology in microglia.
Microglia are involved in many neurodegenerative disorders, such as Parkinson's and Alzheimer's diseases, and are massively activated during brain injuries such as stroke [3,15,16]. There are many known hurdles to overcome in the study of microglia for drug discovery strategies [17][18][19]. Robust and scalable cellular models have to be built that help understand the biology of microglia for drug development. To the best of our knowledge, there is no theoretical work and no successful experimental study of simultaneous activation of several ionotropic receptors in human microglia. There is substantial evidence that both P2X 4 and P2X 7 receptors are expressed in microglial cells [7,20,21]. However, human microglia experiments failed to isolate fucntional P2X 4 receptors [22], therefore a modelling approach will deepen our understanding of the contribution of each of these receptors to Ca 2+ and Na + dynamics within human microglia. Moreover, our model predictions for human microglia and the comparison of the kink-like behaviour of their responses confirm that intricate intracellular Ca 2+ transients arise from a complex cooperative activation of both P2XRs and other ionic extruders and pumps. It is anticipated that the future in situ and in vivo experiments relative to human microglia research will be able to validate the theoretical findings supported by faithful biological inferences made in this paper. There is a high demand to develop reliable cellular human models where microglia can be studied with the aim of furthering our understanding of neuroglia interactions.

Methods
The high-level view of the biophysical model is shown in Fig 1, which includes P2XRs, the Na + /Ca 2+ exchanger (NCX) and the plasma membrane Ca 2+ ATPase (PMCA) along with Ca 2+ and Na + leakage channels. Microglia ionic homeostasis is both dynamic and complex,  2+ . Cellular electrophysiology is included in our model by involving the transmembrane currents that underpin P2X-mediated Ca 2 + signalling and microglia homeostasis. The P2X currents for each receptor are modelled as two individual Na + and Ca 2+ conductance-based channels (I P2X Na and I P2X Ca ), and both Na + and Ca 2+ leak channels (I NaLeak and I CaLeak ) are also included in the proposed model; the dynamic behaviour of the reversal potentials (E Na and E Ca ) are also taken into account. The current through the NCX channel is divided into a sodium current (I NaNCX ) and a calcium current (I CaNCX ).
https://doi.org/10.1371/journal.pcbi.1009520.g001 involving a myriad of receptors, effectors, proteins, channels and processes. The proposed model for intracellular Ca 2+ and Na + dynamics mediated by P2XRs is based on the mass action kinetics formulated by four state variables to capture the complex mechanics of both P2X 4 and P2X 7 receptors when activated by ATP.
Increases in extracellular ATP levels (e.g., from damaged neurons) are detected by microglia and they extend their processes towards the elevated ATP levels through chemotaxis principles [7]. Fig 1A summarises the key intracellular signalling transduction pathways for this purpose. The two P2XRs bind ATP and change their conformation, allowing a rapid rise in intracellular [Ca 2+ ] through Ca 2+ influx. Published experimental data for hP2XRs obtained from both native preparations and expression systems were employed throughout the model development [22,23]. While no data for hP2X 4 R currently exists for microglia, this paper uses functional hP2X 4 data recorded in expression systems [23]. Given the overexpression of receptor protein and the magnitude of expression system currents, the data has been scaled to better align with native microglial P2X currents [24] (in the order of pA).

P2X kinetic models
P2XRs work as cation channels and open after ATP binding. They allow rapid entrance of Na + and Ca 2+ ions into the cell with an efflux of potassium (K + ): note that it has recently been proposed that K + efflux is primarily driven by the TWIK2 K + channel [25] but this mechanism needs further investigation, and therefore K + dynamics are not considered in this article. Both P2X 4 and P2X 7 receptors in microglia share similar trimeric configurations but display different biophysical properties. P2X 7 receptors in the CNS can activate the proliferation of microglia, modulate cell phagocytosis and induce TNF-α production [7]. The P2X 4 receptor, which is stimulated at micromolar ranges of ATP, induces an instantaneous peak current and desensitises within seconds [23,24]. In contrast, millimolar ATP triggers P2X 7 receptors and maintains ionic currents during the ATP application [22].
Biophysically different Markov models have been widely used to describe the binding sites of P2XRs in whole-cell current configurations with eight or more state variables. In this work a new mathematical model was developed that avoids the complex Markov approach [26][27][28]. It is worth noting that eight-state Markov models [29,30] were not successful in providing a good fit to the human data in [22,23] (see S2 Text). Therefore, we developed a minimal model for P2XRs as follows.
Fig 2 shows a biochemical kinetic reaction network that underpins both the hP2X 7 and hP2X 4 receptors comprising of only four state variables (where C, S, D and O states represent closed, sensitised, desensitised and opened respectively). By introducing variable exponential kinetic rates into this four-variable model (inpisred by other similar biophysical models [31]) in a similar way to the approach of Hodgkin and Huxley (HH) [32], the model becomes a less complex realisation as compared to the existing P2X Markov models, which is capable of capturing the major gating properties of microglial P2XRs. Note that the HH model benefits from voltage-gated exponential rates while the model herein makes use of agonist-gated exponential rates. The gating properties of the P2X channels are broken into four individual stages, including: activation (S ! O) and sensitisation (C ! S and D ! S), desensitisation (S ! D and then D ! C) and deactivation (from O to S and then C, or S to D and then to C). Activation is a quick process, which happens after the receptor is sensitized by an agonist, during which the channel opens resulting in increasing inwardly cationic currents following ATP exposure which slowly tends towards a plateau. This leveling off of the current is reffered to as desensitisation and following this the current amplitude rapidly goes to zero when the ATP becomes unbound. P2X 4 and P2X 7 receptors differ in terms of the kinetic machinery actively involved in their sensitivity and function.
In the model, ATP application is potentiated though two forward rate constants of α 1 and α 2 . All backward rates were chosen as variable functions. Experimental data show that desensitisation of P2XRs depends on the level of ATP exposure [23,29,33,34]. Therefore, the gating transition between S⇄D has an exponential dependency on the agonist. This prevents the desensitisation pathway from disappearing in the model during agonist removal since the exponential terms became unity upon ATP removal. Note that both exponential rates are present in both sides of the transitions guaranteeing the amount of reduction in S state in the forward reaction returns to D state in a proportionally balanced fashion. The model conserves all the species because it is a closed system; therefore the the sum of all state variables at any time is equal to unity and the system state variables return to their initial resting state after stimulation has ceased. An important function of a biophysical model is that the state of the system should go back to its resting state, because the model can be activated again without affecting the model predictions, for example, when this model is simulated using a repeated application of ATP. In order to ensure this property, the transition between D and C states through a negative exponential rate function makes an extra path that enables D to be depleted to C when O becomes very small (note that the path is activated because the exponential function goes to its high value when O significantly decreases).
The microglial hP2X 7 receptor currents [22] have two significant phases: A to B (activation/ sensitisation/desensitisation) and B to C (deactivation). The exponential backward rates in Fig 2 between O/S and S/C allows this complex behaviour to be captured. The function φ for the hP2X 7 model is defined as φ(A) = A. The exponential backward rates in Fig 2 between O/S and S/C makes this complex behaviour possible to be captured in such a simple manner. In fact, when agonist is high the exponential terms are low and let the model pass the flow in the forward direction (C/S/O). By contrast, these two terms switch to their maximal values (viz. one) upon agonist removal, so letting the state O deplete (go to zero) in order to capture the deactivation phase (B to C) in Fig 3A. Strictly speaking, such intricate dynamics can be captured by simpler models when reaction rates are desiged such that they dynamically change to avail of existing paths of the model. They are responsible for reproducing the main properties of the experimental data for whole-cell current, such as persistent sensitised state and monophasic/biphasic dynamics.
In contrast, hP2X 4 receptors mediate a rapid current that is sustainably greater than the peak amplitude of P2X 7 receptor as shown in Fig 3B, and their desensitisation is much faster https://doi.org/10.1371/journal.pcbi.1009520.g002 than P2X 7 R, which is also validated by the model predictions. As seen in Fig 3B, hP2X 4 receptor has three significant phases at macroscopic level: A to B (activation/sensitisation), B to C (desensitisation), and C to D (deactivation). In contrast to microglial rP2X4 [24], the hP2X 4 current takes some time to vanish after ATP removal (C to D). This additional phase requires minor changes to the P2X model of Fig 2, where φðAÞ ¼ e z 1 A is chosen. The transition between S and D states are gated by an exponential form of the agonist in Fig 2. It allows the model does not need additional state variables and therefore avoids the necessity for a more complex model. When ATP is zero, the transition between C and S is eliminated while φ(A) becomes unity. In fact, the bidirectional path between S/D and S/O makes active routes for the current to vanish, namely, the phase of C to D in Fig 3B is  To derive the mathematical equations representing the kinetic P2X model, the topology in Fig 2 is now translated into a system of time-dependent coupled ordinary differential equations as expressed in Eqs (1-4) by using the law of mass action: https://doi.org/10.1371/journal.pcbi.1009520.g003 Each state C, S, D and O corresponds to the fraction of a specific P2X receptor that are in the Closed, Sensitised, Desensitised and Opened state at any given time. The above equations are used for both P2X models (P2X 4 and P2X 7 ), where Z is used to denote the subtype receptor number (4 or 7). The fitting of the model parameters is discussed in S1 Text. Values and units of rate constants and coefficients in Eqs (1)(2)(3)(4) are provided in the model fitting and validation section.

hP2X Z current model
Having the system of Eqs (1) to (4), we can now derive the current equations for the P2XRs (4 and 7). Currents through the receptors are proportional to the state O Z , where Z is the subchannel number of the receptor, and can be expressed as terms of the familiar conductance channel model. These currents are obtained by a product of the open state, the maximum conductance g through the channels, and the difference between the membrane potential (V m ) and the reversal potential of the cell (E). The Ca 2+ and Na + currents for the receptors can be modelled as two separate conductance channels using Eqs (5) and (6), as illustrated in Fig 1B, where reversal potentials are expressed in Eqs (8) and (9). With regards to K + efflux, the mechanisms underpinning this ionic current are poorly understood and are therefore not considered in this paper. However, the influx of Na + will affect the NCX and will consequently have a downstream effect on Ca 2+ . Therefore, the contribution to the Na + concentration in the cytoplasm by the P2XRs needs to be considered for a more complete understanding of the observed Na + and Ca 2+ dynamics. This approach, where a conductance type model is used to characterise the P2XRs as individually ATP-gated Ca 2+ and Na + channels, allows the segregation of the Ca 2+ or Na + channel currents while maintaining a similar overall time-dependant response. While we recognise that this is an assumption, the segregation of these two ions facilitate the calculation of both Ca 2 + or Na + ionic concentrations in the cytoplasm. Ca 2+ or Na + channel conductance in the proposed model is set to define channel stoichiometry. The total current is calculated by the summation of the separate channel contributions, as shown in Eq (7).
The above equations are used for both P2X models (P2X 4 and P2X 7 ), where Z is used to denote the subtype receptor number (4 or 7). The ionic concentrations in the extracellular space are assumed constant but are dynamic in the cytoplasm due to the influx/efflux of ions giving rise to a dynamic transmembrane concentration gradient. Accordingly, the reversal potentials for Na + and Ca 2+ in Eqs (5) and (6) can be expressed in terms of the Nernst equation as where x and i are are the extracellular and intracellular space respectively, Z Na þ and Z Ca 2þ are the valencies of Na + and Ca 2+ respectiveley, F is Faradays constant, R is Boltzmann constant and T is temperature in kelvins.

Sodium calcium exchanger (NCX)
The NCX is an antiporter that electrochemically exchanges three Na + ions for one Ca 2+ ion across the cell membrane [35]. This exchanger works in two different operational modes (forward or reversed) depending on V m and the transmembrane Na + and Ca 2+ gradients. The Na + and Ca 2+ current densities of the exchanger are defined by where � J NCX and γ are the NCX conductance and partition number, respectively [35]. The current density J NaNCX has been obtained through a membrane-coupled protein that assumes both Na + and Ca 2+ are abundant in the extracellular space. The current density therefore depends on the intracellular and extracellular concentrations of ions and V m . The exchanger has a low affinity and binds moderately to Ca 2+ . It transmits both Na + and Ca 2+ ions quickly because of its high capacity and consequently the NCX is mainly responsible for removing Ca 2 + that flows into the cytosol.

Plasma membrane Ca 2+ ATPase (PMCA)
The PMCA pump is a high-affinity protein that preserves low cytosolic Ca 2+ profiles by consuming ATP macromolecules and extrudes Ca 2+ from inside the microglia against its ionic gradient. The Ca 2+ current density through PMCA is modelled using Michaelis-Menten kinetics [36], and is formulated by where � J PMCA and K PMCA are the maximum current density and pump affinity respectively. n is the order of Hill function, which is set to one for this study. This pump removes Ca 2+ and is also used to maintain Ca 2+ at low levels in our model.

Na + and Ca 2+ leakage channels
In our model, two separate membrane bound leakage channels exist for Na + and Ca 2+ and are modelled as two passive electrochemical gradients by assuming a Nernst-like equation given by where g y and E y are a conductance and a reversal potential for ion y, and V m is the membrane potential. The physiological conditions of these channels change dynamically with both membrane and reversal potentials.

Na + and Ca 2+ concentration dynamics
Na + and Ca 2+ signalling is formulated as a fluid compartment model, which is assumed to be well mixed. The changes in ½Na þ i � and ½Ca 2þ i � are modelled as where S m , S p , Z, F and V p are the microglial membrane surface area, microglial process surface area, the valency of ions, Faradays constant and the intracellular process volume respectively: the model assumes that ½Na þ i � and ½Ca 2þ i � reside in the total cytosolic volume of microglial processes. The conversion of the P2XR's whole-cell patch-clamp currents into current densities is performed by dividing over S m due to other current components already being expressed as current densities.

Membrane potential
Taking the same approach as is used for neuronal cells [32,37], V m can be related to ionic flow across the membrane as where C m is the membrane capacitance. The total current flow through P2X channels are added to the transmembrane currents across NCX exchanger, PMCA pump and leak channels.
The terms expressed as current densities in Eq (17) are multiplied by S p .

Model fitting and validation
In this section, we fit (see S1 Text) and validate our P2X models using available experimental data [22,23]. As illustrated in Fig 3, the kinetic P2X model shown in Fig 2 is able to reproduce all the major gating properties of hP2X 7 and hP2X 4 receptors for the whole-cell current experiments, i.e. activation, deactivation, sensitisation and desensitisation. The fitted hP2X 7 model is subject to 5mM ATP stimulation [22], while 0.1mM ATP had been used for the hP2X 4 receptor [23]. In contrast to hP2X 4 current profile that peaks at -372pA, hP2X 7 responds much slower reaching a peak of -145pA. Clearly, both models closely approximate the experimental data. The right hand side of Fig 3 illustrates transient behaviour of the whole-cell current state variables. As seen, sensitisation and desensitisation significantly differ in both models. Interestingly, there are two similar transients in state S of Fig 3A that align with the the existing two phases of activation and deactivation. The first transient relates to the activation phase from point A to B (Fig 3A) where ATP is present, and the second transient is created upon ATP removal since the state O must be depleted from the only available path (O to S, followed by S to C or S to D and then finally to C in Fig 2A). The summation of transient responses of all the state variables are equal to unity at all times and all state variables return to their resting state post stimulation. There are four classes of parameters in the entire model formulated by Eqs (1) to (17), including, (a) whole-cell patch-clamp current parameters of α and β, (b) estimated constant parameters such as z, (c) morphological parameters such as C m and S p , and (d) physically measured parameters such as K PMCA and γ used from the literature. Conductances of both P2X models were fitted to clamped electrophysiological voltage data [22,23] and then reduced by a factor of 25 (similar to [30]) such that the difference between the resting [Ca 2+ ] and the peak [Ca 2+ ] (as can been seen later inresults section) lies between 50nM and 55nM, in agreement with the general P2X-mediated ½Ca 2þ i � dynamics found in the literature [38]. PMCA, NCX and leak channel parameters were determined through numerical simulation to fit experimental data [39] such that ½Na þ i � and ½Ca 2þ i � are maintained at their resting concentrations in steady state. This is a paramount property of the microglial cells where the model takes up the role of maintaining ½Ca 2þ i � in resting microglia because microglia must be sufficiently sensitive to their microenvironment to generate a quick response to injury. By letting the optimiser calibrate the value of z parameters in Fig 2, it was unable to give a satisfactory fit quality. Hence, z was estimated manually by setting different values for z and letting the optimiser find all remaining parameters. Therefore, trial and error was used where values of z 1 and z 2 greater than 1 was chosen, and these values were then used in the optimiser. This process was repeated until the best fit for all remianing parameters was found (see Table 1 for final values of z 1 and z 2 ). S5 Text gives an example for choosing a typical value of z 1 in that the model cannot capture the experimental data well.
Morphological changes occur in microglia and while the proposed model does not take this into account, morphological data [40] was used to estimate the surface area and volume of activated microglia by assuming three cylindrical processes and the cell body is represented by a cube: while microglia can dynamically lose their spherical somatic phenotype upon activation, a cube is closer to an image of particular microglia morphology as reported in [40]. Using these assumptions, the process surface area and volume can be calculated and used as parameters in Eqs (15) to (17). Values, units and descriptions of the parameters used in our model are presented in Table 1.

Results
The model developed in the previous section is now studied by numerical simulation to provide a deeper insight into intracellular P2X-triggered Na + and Ca 2+ dynamics in microglial cells. The model was implemented in MATLAB Release 2021a. We used ode23s routine to numerically integrate the non-linear, stiff model equations. The integration was carried out using the default MATLAB ODE solver timestep because it takes advantage of dynamic step sizes for implicit numerical integration. For integration stability during the curve fitting and numerical simulation, the calculation of the right-hand-side terms of the ODEs in Eqs (1) through (17) and the Jacobian matrix of the system equations were implemented for the model. This section investigates four case studies, including, whole-cell current due to P2XRs, Ca 2+ , Na + and V m dynamics resulting from the P2X-provoked currents.

P2X-mediated whole-cell currents
Microglial hP2X 4 and hP2X 7 receptors are activated in micromolar and millimolar concentrations of the applied ATP, respectively [24,38]. We now consider model predictions for both receptors. The simulated whole-cell current responses of hP2X 7 and hP2X 4 receptors are shown in Figs 4-6 for different applications of ATP stimulus level and duration.
All the results for hP2X 7 receptor invoked by micromolar or millimolar ATP (see Figs 4 and 5 respectively) display two biphasically distinguishable phases as observed in [22], namely,  [22]. Note that: it is inferred that ATP levels of 1-2mM corresponds to 300μM of Benzoylbenzoyl-ATP (Bz-ATP) [41]. This is also discussed further in the next section. The predictions in Fig 5A and 5B show that for short ATP application, the current remains in ascendancy while for long ATP stimulus the current plateaus and thereafter slowly declines. In contrast to rat microglial data [42], the current amplitude appears to saturate with larger ATP levels, which suggests that desensitisation becomes more pronounced for higher levels of ATP. However, this behaviour agrees with what has been reported for mouse P2X 7 microglia receptors [43].
The predictions in Fig 6A show that the amplitude and rate of increase of current is largely unaffected by the magnitude of the ATP stimulus for hP2X 4 R. Note that this receptor appears to rapidly desensitise until the ATP stimulus is removed, thereafter the current falls off rapidly to baseline. In contrast, Fig 6B shows that for long ATP durations the current reaches a steady state after a period of desensitisation and these levels are ATP dependant. Only when the ATP  stimulus is removed does the current return to its baseline level. This dual component indicates a biophsyical fingerprint that has been reported through in vitro electrophysiology [44]. For example, human monocyte-derived macrophage research shows that P2X 4 Rs exhibit fast activation current followed by biphasic desensitisation in [34], particularly for high levels of ATP. Pertinent to dynamic microglia, this biphasic fucntionality of P2X 4 Rs may provide a regulatory pathway to control motility as has been suggested in expression systems [27].
The hP2X 7 and hP2X 4 receptor models can now be integrated yielding valuable data regarding P2X-mediated Ca 2+ and Na + signalling in human microglia. These data could be used to guide future experimental design.

P2X-mediated calcium dynamics
The interaction between the P2XRs, NCX, PMCA and leak channels dictates the dynamics of Ca 2+ and Na + concentrations and also V m , as described by Eqs (15) to (17). The model predicts the total Ca 2+ influx/efflux via the receptors, leak channels, extruders and pumps. Fig 7 shows individual Ca 2+ transients for each receptor separately assuming the same ATP exposure level and duration. For Fig 7A, the ½Ca 2þ i � reaches a plateau for higher levels of ATP and this correlates with the total current carried by this receptor (e.g., at 3mM), as shown in Fig 5B. Fig 7B  agrees with experimental intracellular Ca 2+ measurements [22] for ATP = 1mM, as discussed earlier, in which there exists a local maximum before the Ca 2+ curve goes to its plateau prior to ATP removal. Fig 7B shows a very transient-like ½Ca 2þ i � profile where hP2X 4 R preserves ½Ca 2þ i � for higher levels of ATP which is consistent with its equivalent current profile in Fig 6.  Fig 8 shows the absolute currents carried by the P2X 7 (P2X 4 switched off), NCX, PMCA and leak channels for ATP levels of 0.1mM, 1mM and 3mM. Note that prior to the onset of the ATP stimulus Ca 2+ efflux by the PMCA is negated by the NCX operating in reverse mode and calcium leak across the membrane via the Ca 2+ leak channel. As the ATP stimulus is applied, the P2X 7 receptor is activated and both Ca 2+ and Na + current flows into the cytoplasm. However, V m and ½Ca 2þ i � are the primary drivers for the NCX ensuring that it operates in reverse mode, as can be seen in Fig 8 for higher levels of ATP (1mM and 3mM). Also, after the initial disturbance, both the Ca 2+ and Na + currents stabilise such that the net Ca 2+ or Na + current Reversal potentials and V m for P2X 7 R are shown in Fig 9 when P2X 4 R is switched off. Note that the reversal potential for Na + (Fig 9) changes much less significantly in comparison with ½Ca 2þ i �, because the baseline level of ½Na þ i � is much higher than that of ½Ca 2þ i �. Additionally, experimental data [42] indicates that rectification of V m occurs following low levels of ATP  Fig 9C) shows a similar trend to this data.
Note that the mode of operation for the NCX is, to a good approximation, unaffected by influx/efflux of Na + currents since the baseline ½Na þ i � in the cytoplasm is high (8mM). In contrast, the mode of operation of the NCX is much more sensitive to the ½Ca 2þ i � in the cytoplasm as this quantity changes much more significantly by the ATP stimulus because the baseline ½Ca 2þ i � in the cytoplasm is low (45nM); V m also drives the NCX. Fig 10 shows the absolute currents carried by the P2X 4 , NCX, PMCA and leak channels for ATP levels of 0.1mM, 1mM and 3mM. Prior to the onset of the ATP stimulus the Na + leak current is negated by NCX extruded Na + efflux and the Ca 2+ efflux by the PMCA is negated by the NCX operating in reverse mode and calcium leak across the membrane via the Ca 2+ leak channel. As the ATP stimulus is applied, the P2X 4 receptor is activated and both Ca 2+ and Na + current flows into the cytoplasm. In contrast to the P2X 7 current (see Fig 5), the current through the P2X 4 (see Fig 7) is more transient in nature, reflecting the behaviour of the biological P2X 4 receptor. The sharp rise in Na + influx by the P2X 4 results in an associated rise in the Na + leakage current and V m , with the NCX extruder driven deeper into reverse mode. Also, the sharp rise in Ca 2+ influx resulting from P2X 4 activation causes additional efflux of Ca 2+ via the PMCA but this is somewhat negated by the influx of Ca 2+ by the NCX, where the latter is driven by hyperpolarisation of V m and fluctuations in the ½Ca 2þ i � in the cytoplasm. As in the case of the P2X 7 , the reversal potential for Na + (Fig 11) is virtually unaffected due to the high ½Na þ i � baseline in the cytoplasm but the reversal potential for Ca 2+ (Fig 11) changes much more significantly because of the low ½Ca 2þ i � baseline in the cytoplasm. A holistic view of ½Ca 2þ i � transients by concurrent stimulation of both hP2X 4 and hP2X 7 receptors is given in Fig 12. A significantly distinguished pattern is observed in Fig 12 for ATP = 1mM and ATP = 3mM. As seen in the results of Fig 8, both the Ca 2+ and Na + currents continue towards saturation while the ATP stimulus is sustained. In contrast, for the P2X 4 (Fig 10) both the Ca 2+ and Na + current behave more like a transient (effecting V m ) and decay over time irrespective of the level or duration of the ATP stimulus. It is therefore concluded that the kink in the curves at ATP levels of 1mM and 3mM in Fig 12 are because of the fast fall off in ½Ca 2þ i � due to the decrease in the Ca 2+ current through the P2X 4 and the relatively slow rise in ½Ca 2þ i � due to the increasing Ca 2+ current through the P2X 7 . When these two ½Ca 2þ i � cross over (at a comparable concentration) a kink is formed (momentary rise in ½Ca 2þ i � followed by a slow fall off in ½Ca 2þ i �). Essentially, both receptors have direct influence on the total P2X-mediated ½Ca 2þ i � dynamics where P2X 4 receptor dominates the ½Ca 2þ i � at the onset of the ATP stimulus but thereafter the P2X 7 receptor maintains a high level of ½Ca 2þ i � in the cytoplasm until the removal of the ATP stimulus. A wide spectrum of experimental data confirms this theoretical finding where several P2 (i.e. P2X and P2Y variants) receptors are simultaneously triggered by extracellular nucleotides like ATP. Activation of P2X 4 R and P2X 7 R results in intricate ½Ca 2þ i � kinks for 1mM and 3mM of ATP in rat microglia [38]. Other work on microglia and astrocytes reveal kink-style behaviour of ½Ca 2þ i � [16,[45][46][47] and such dynamics are also observed in other animal and human cell lines including those experimentally reported in [48][49][50][51][52][53]. Note that no kink appears for ATP < 1mM because the initial rise in current is dominated by current flow through the P2X 7 . Fig 13 shows the ½Na þ i � and V m dynamics. As is evident, both ½Na þ i � and V m graphs have similar shapes in contrast to ½Ca 2þ i � counterparts, as their dynamics are dictated by the action of the P2XRs. It is because NCX channels mainly control the flow of Na + cations across the membrane as microglia are reported not to express Na + /K + -ATPases [7]. The behaviour of V m indicate that microglia depolarise rapidly, which is also consistent with experimentally reported data for rate microglia [42]. Emerging evidence shows that the dynamic membrane potential, V m , plays key roles in regulation/upregulation of biological processes such as cell migration and proliferation [54]. Therefore, the model may be made more biologically plausible by considering voltage-gated Na + channels and ½K þ i � dynamics. It is worth noting that varying the morphological parameters of the model such as S p and V p only alters the predicted ½Ca 2þ i � amplitudes slightly and so does not significantly influence the overall response of the Ca 2+ dynamics (see S3 Text). In contrast to the current profile outlined in Fig 8 for the hP2X 7 only and in Fig 10 for the hP2X 4 only, the marked difference when both receptors are activated simultaneously (as seen in Fig 14), is the domination of the hP2X 4 receptor current at high levels of ATP, while the hP2X 7 stays constant at higher ATP levels and over a longer duration. In addition the kink in the concentration profile for Ca 2 directly correlates with the initial fall off in the Ca 2+ current, carried by the PMCA, due to the decaying Ca 2+ influx by hP2X 4 and the increasing Ca 2+ influx by hP2X 7 (see curves for 1mM and 3mM ATP levels).
To summarise, it is inferred that the very transient behaviour of the kinks in Fig 12 is indicative of cooperative dynamics between activation of several P2XRs, PMCA and NCX, where NCX mainly derives the ½Na þ i � extrusion. As such, V m is similar to Na + dynamic in Fig 13  because of a high baseline of ½Na þ i �, and leak channels keep microglia sensitive enough to changes in their microenvironment.

Sensitivity analysis of the mathematical model
As demonstrated throughout this paper P2X-mediated Ca 2+ dynamics are primarily controlled by ATP-gated P2XRs, NCX, PMCA and leak channels. This section carries out sensitivity analysis (SA) of the model parameters to elucidate the robustness of the model responses [55]. It is also a powerful means whereby the parameters or components that substantially contribute to the model predictions can be determined. The sensitivity analysis results can provide a detailed insight into the model parameters for future model refinements when new experimental data becomes available in order to calibrate the model at different levels of agonist. Sensitivity can also be viewed as the inverse of robustness where sensitivity of parameters give rise to a quantitative estimate in deviations of model outcomes arising from perturbation of model parameters. Since the behaviour of high-dimensional biochemical reaction networks is often dominated by a small number of parameters, SA helps identify these parameters that can be investigated in further studies.
In this section, two different methods of SA ar are employed to obtain additional information regarding the P2X model parameters: these methods are local and global sensitivity analysis (LSA and GSA). In LSA, a single parameter (where others are kept constant) is varied and the effect of this perturbation is considered. This process is repeated for all parameters. In contrast to LSA, GSA techniques employ simultaneous pertubations to model parameters, and changes in the model output are monitored. GSA is capable of detecting relationships between a set of parameters. Logarithmic sensitivity analysis (LSA) [56], was initially used to investigate the influence of the model parameters on the model output. LSA is dimensionless and therefore allows comparison of the results conveniently. If y is a single response of a system of ODEs, LSA of this variable is defined as where p i denotes the i'th parameter. Partial derivative of y in Eq (18) with respect to p i can be approximated by using forward finite difference (FD) [55] as follows where Δp i is a notation for a small change of the parameter p i . As a general rule of understanding sensitivity analysis robustness, if the results do not change significantly after the fitted model parameters are perturbed then the sensitivity analysis is said to be robust [55]. For this study, the sensitivity analysis was carried out on the cytosolic Ca 2+ concentration with respect to the maximum NCX current density ( � J NCX ), the maximum PMCA current density ( � J PMCA ), and the parameter sets of the P2X models by Several techniques exist to carry out GSA. Sobol's method is used in this section, which is a variance-based method [57]. The method does not make any assumption for input and output relationships of a model. The variance of the model output is decomposed into a combination of variances. A Joint Probability Function (PDF) is employed to make the Sobol's terms. A MATLAB implementation of this method developed by [58] was used to perform GSA and study the effect of simultaneous changes of model parameters on cytosolic Ca 2+ for ATP = 1mM. Fig 16 shows the computed Sobol's total indices for the peak duration. A peak duration is defined as the time interval at which the transient values of the model output is greater than or equal to the mean of the output throughout the simulation time. According to Fig 16 for total Sobol indices, the five most dominant parameters are α 3x4 , α 2x4, β 2x4 , β 2x7 and β 3x4 in decreasing order. β 2x7 and g x7 have more influence on the output among all other parameters of the P2X 7 model. The effect of both P2X conducntances is nearly close, while there is a significant variation for most of the corresponding parameters of both P2X models. β 2x7 and β 2x4 have the biggest impact on the output among all other backward rate coefficients for each individual model. Interestingly, the global sensitivity of g x7 in Fig 16 is significantly less influential than its local sensitivity illustrated in Fig 15. These LSA and GSA results show that the parameters of the P2X 4 R model are mostly dominant in the model sensitivity as a whole.

Discussion
Neurotransmitters can modulate intracellular Ca 2+ and Na + via activation of microglial receptors, which regulate their homeostatic mechanisms [7]. Microglia express P2X 4 and P2X 7 purinergic receptors which impact local Ca 2+ homeostasis. Mathematical models that result in a formal description of agonist binding to purinergic receptors with receptor activation, provide unique tools to unravel the underlying regulatory dynamics of intracellular ionic concentrations. In this research a minimal computational model for Ca 2+ , Na + and V m dynamics of microglia was developed. The binding of ATP to P2XR on microglial processes leads to receptor activation, with the associated influx of Ca 2+ and Na + . This results in a rise mainly in ½Ca 2þ i � and V m hyperpolarising, leading to activation of the PMCA, NCX and leak channels: it is worthwhile noting that this model allows a detailed investigation into the interplay between the P2XRs, pumps, extruders, membrane voltage and leak channels. The model predictions show that activation and deactivation of P2X 4 and P2X 7 receptors in microglia generate monophasic and biphasic shapes in conformance to microglia-specific P2X data or other cell-specific data at lower and higher ATP concentrations. Particularly, for the P2X 7 receptor, there are two different phases of current initialisation and current facilitation. Simulation observations show that monophasic and biphasic transients are also present in single cell Ca 2+ and Na + profiles. The model indicates that P2X 7 receptor maintains significant currents and therefore effect ionic concentrations much more than P2X 4 at higher levels of ATP. These data have implications for brain pathology, where elevated levels of ATP are evident [59,60]. The relative contributions of the P2XRs studied here may provide insight into therapeutic strategies targeting either receptor. For example, P2X 7 antagonism has been suggested as a therapeutic intervention for a range of CNS disorders [61][62][63]; here the modelling data present a read-out of this pharmacological intervention in relation to ionic homeostasis and microglia function. In addition, there are implications from the model with respect to microglial function where their levels of the P2XRs are known to be altered (e.g. Parkinson's disease; [64]).
The model demonstrates that Ca 2+ dynamics not only rely on different gating property of P2XRs as a function of ATP pulse protocol but also on the opposing Ca 2+ currents carried by PMCA and NCX: the latter being continually in the forward mode due to the hyperpolarisation of V m and ½Ca 2þ i � fluctuations in the cytoplasm. At ATP levels of 1mM and 3mM, unexpected kinks are observed in the concentration of Ca 2+ regulated by ½Na þ i � and V m . This kink is a direct consequence of the fast fall off in ½Ca 2þ i � due to the continual decrease in the Ca 2+ influx, via the P2X 4 competing with the slow rise in ½Ca 2þ i � due to the increasing Ca 2+ influx via the P2X 7 . Hence there exists a cross over in the ½Ca 2þ i � and a kink appears followed by a gradual decrease in ½Ca 2þ i �: the P2X 4 receptor dominates the ½Ca 2þ i � at the onset of the ATP stimulus and thereafter the P2X 7 receptor maintains a high level of ½Ca 2þ i � in the cytoplasm until the removal of the stimulus. The abstraction level of the proposed model facilitates the reproduction of meaningful responses arising from the interplay between the P2XRs, PMCA, NCX and leak channels.
The current experimental data shows a total cellular current where the contribution to the total current comes from Ca 2+ and Na + influx, and K + efflux. With regards to the latter, the mechanisms underpinning K + efflux are poorly understood and are therefore not considered in this paper. However, the influx of Na + will affect the NCX and therefore have a downstream effect on Ca 2+ . Because of this, the contribution to the Na + concentration in the cytoplasm by the P2XRs was considered for a better understanding of the observed Na + and Ca 2+ dynamics. To the best of our knowledge, there is no theoretical or experimental study of activation of multiple ionotropic receptors in human microglia. Cultured human microglia experiments failed to isolate P2X 4 currents [22] nor did they provide data on the interaction between different P2XRs. However, the model presented in this paper does provide insight into both of these receptors and how they collectively contribute in terms of Ca 2+ and Na + dynamics. Moreover, the simulation results in this paper used human data in so far as possible and can inform experimentalists as to the role of Na + influx and how this can affect downstream Ca 2+ dynamics.
The proposed computational framework provides a strong foundation that can be extended towards a more complete understanding of microglial function, particularly, in situ and in vivo. The model can also be refined as new experimental data becomes available and to this end the expansion in human microglia research through the use of induced pluripotent stem cells is of interest [65]. Additionally, the function of the P2Y 12 receptor will be the focus of future work, as these are highly expressed in activated microglia [66]. It has been shown that P2Y 12 actively participate in triggering the PI3K pathway through Ca 2+ uptake from the extracellular space [8,67]. Other future work on extending the model proposed here is the inclusion of distinct types of K + , and voltage-gated Na + channels, which have rapid kinetics and depend on V m . The biophysical model developed in this paper will improve our quantitative understanding of P2X-mediated microglial physiology.