Deciphering the regulation of P2X4 receptor channel gating by ivermectin using Markov models

The P2X4 receptor (P2X4R) is a member of a family of purinergic channels activated by extracellular ATP through three orthosteric binding sites and allosterically regulated by ivermectin (IVM), a broad-spectrum antiparasitic agent. Treatment with IVM increases the efficacy of ATP to activate P2X4R, slows both receptor desensitization during sustained ATP application and receptor deactivation after ATP washout, and makes the receptor pore permeable to NMDG+, a large organic cation. Previously, we developed a Markov model based on the presence of one IVM binding site, which described some effects of IVM on rat P2X4R. Here we present two novel models, both with three IVM binding sites. The simpler one-layer model can reproduce many of the observed time series of evoked currents, but does not capture well the short time scales of activation, desensitization, and deactivation. A more complex two-layer model can reproduce the transient changes in desensitization observed upon IVM application, the significant increase in ATP-induced current amplitudes at low IVM concentrations, and the modest increase in the unitary conductance. In addition, the two-layer model suggests that this receptor can exist in a deeply inactivated state, not responsive to ATP, and that its desensitization rate can be altered by each of the three IVM binding sites. In summary, this study provides a detailed analysis of P2X4R kinetics and elucidates the orthosteric and allosteric mechanisms regulating its channel gating.


Introduction
Purinergic P2X receptors (P2XRs) are a family of ligand-gated non-selective cation channels that are activated by extracellular adenosine 5'-triphosphate (ATP). In mammals, there are seven distinct subunits of this family of proteins, labeled P2X1-7. Each subunit contains intracellular N-and C-termini connected to the first and second transmembrane (TM) domains, respectively, followed by a large extracellular loop commonly referred to as the ectodomain. It is also well established that P2X subunits aggregate to form functional trimers [1][2][3][4][5]; receptors may be composed of either one type of subunit (homotrimer) or a mixture of more than one type of subunit (heterotrimer) [6]. When coordinated in a trimer, the interfaces between adjacent ectodomains form three binding pockets for ATP [7]. These ectodomains also form fenestrations which are lined by negatively charged amino acids that attract cations, and the cation selectivity of these channels is determined by the selectivity filter localized in the TM2 domain [8].
The binding of two or three ATP molecules to the extracellular binding sites induces conformational changes in the ectodomain and subsequently the TM domains, causing the channel opening. The gating of P2XR by ATP and other orthosteric agonists can be broken down into three distinguishable phases, activation, desensitization, and deactivation, defined by their ionic current kinetics in whole-cell recordings [9,10]. Activation is a rapid phase of channel opening that corresponds to increasing inward current subsequent to agonist application. This is usually followed by the desensitization phase, a decay of current amplitude in the maintained presence of an agonist, with an onset that is slower than that of activation. After agonist removal from the medium, a relatively rapid decrease in current amplitude, referred to as the deactivation phase, is observed. Receptors differ in both their sensitivity to agonists and the kinetics of the phases (with desensitization being transient and controversial in P2X7Rs) [9,11].
It has also been suggested that P2X2R and P2X7R are capable of exhibiting another phase in their gating, termed dilation, when the receptor pore was thought to progressively enlarge during sustained ATP application. Two observations were used as evidence for pore dilation: the ability of these receptors to become permeable to N-methyl-D-glucamine (NMDG + ), a large organic cation (~7.3 Å in mean diameter), and the change in reversal potential (E rev ) during a ramp protocol, when cells were bathed in a medium containing only NMDG + with pipette containing only NaCl [12]. Recent investigations, however, have shown that changes in E rev during prolonged channel activation of P2X2R do not reflect pore dilation, but rather time-dependent alterations in the concentration of intracellular ions, specifically washout of intracellular Na + and gain of NMDG + through the initially opened large pore of P2XR [13,14] and that permeation of NMDG + can also occur without pore dilation [15].
In addition to orthosteric regulation, P2XRs also exhibit allosteric regulation [9]), as is evident from the action of ivermectin (IVM) on P2X4R channels. Extracellularly applied IVM increases current amplitude at low concentrations, and increases the sensitivity of receptors to ATP and partial agonists at higher concentrations. IVM also decreases the extent of desensitization in the continuous presence of agonist and prolongs deactivation of the receptor after the removal of agonists [16][17][18]. Furthermore, P2X4R is not substantially permeable to NMDG + natively [16,19], but displays a shift in E rev in the presence of IVM, suggesting that the channel pore is permeable to NMDG + in IVM-treated cells [20].
The action of IVM on P2X4R gating is also time-dependent; i.e. the cells must be exposed to IVM for at least 30 s, in the absence of ATP, to alter the P2X4R gating (compared to ms for orthosteric activation) [20]. The onset of IVM's potentiating effect on P2X4R current amplitude is faster than the effects of IVM on deactivation kinetics [17,21]. Consequently it has been postulated that the two distinct effects of IVM are due to binding at two distinct sites [17]. Experiments with chimeric receptors containing domains from IVM-sensitive P2X4R and the IVM-insensitive P2X2R have provided evidence that TM domains play a critical role in this allosteric modulation by IVM [18]. The location of the IVM binding site has not yet been addressed in the context of the recent crystal structures of a zfP2X4R [2,3]. However, IVM apparently inserts between pairs of neighboring subunits of the P2X4R channel in the membrane and interferes with the molecular rearrangement in the TM domains involved in channel gating, similarly to glutamate-gated chloride channels crystalized with IVMs [22]. Accordingly, there should be three potential binding sites for IVM in the P2X4R, as there are three clefts between subunits. Such a topography of IVM binding sites provides rationale as to why receptors (not previously stimulated orthosterically) must be exposed to IVM for a prolonged period for it to be effective. Subsequently, we have used the term priming to describe the time and concentration dependence of IVM to occupy its binding sites and the resulting development of its varied allosteric effects.
In recent years, mathematical modeling has begun to shed light on many aspects of P2XRs and to guide experimental designs to arrive at a more complete understanding of channel gating [19,20,[23][24][25]. Biophysically detailed Markov models that describe individual orthosteric binding sites and their allosteric modulation, have been very successful in deciphering the kinetics of P2X homotrimers and succinctly explaining many phenomena [19,20,[24][25][26]. They consider important biophysical details such as the conformational states of individual binding sites and other structural components of the receptor. One of these models for P2X4R is a simple Markov model that takes into account the sequential binding of ATP to its three subunits and assumes that IVM causes receptor sensitization upon the binding of three ATP molecules, that all ATP unbinding rates are decreasing functions of IVM concentration, and that IVM induces a change in ion selectivity caused by the assumed pore dilation [20]. In this model, the three allosteric effects of IVM on P2X4Rs are induced by a single IVM-dependent transition that allowed for generating the shift in E rev during the ramp protocol. However, the published model is unable to account for effects of pre-treatment with IVM before ATP application. The model also predicts a large (> 150%) increase in the unitary (single-channel) conductance of the receptor, in contrast to experimental evidence [17] indicating that there is at most 20% increase in unitary conductance.
To satisfy these constraints, we developed two substantially larger models that not only fit the data more closely in more experimental circumstances but offer better insights into how the kinetics of ATP and IVM sequential binding to P2X4R affect P2X4R activation, desensitization, and deactivation. They also illustrate how changes in ion selectivity of these receptors are manifested, as well as predict the previously unappreciated existence of receptor states (including the deeply inactivated and primed states) that are not directly observable in the experimental current recordings.

IVM reduces P2X4R pore selectivity
When HEK293 cells expressing rat P2X4R were bathed in Ca 2+ -containing medium, where Na + was substituted by NMDG + and a voltage ramp from −80 mV to +80 mV was applied to the cell, E rev was not found to change during sustained applications of 100 μM concentrations of ATP (S1 Fig, left). However, pretreatment with 3 μM IVM for 60 s caused a positive shift in E rev during sustained applications of 100 μM ATP (S1 Fig, right). This is consistent with our earlier work [20] and the finding of others that IVM potentiates ATP-induced responses and increases permeability for NMDG + , but cannot activate P2X4R channels on its own [16]. Because strategies that rely on changes in E rev to provide evidence for large pore formation during sustained stimulation with agonist were questioned [13,15], we examined currents induced by ATP in Ca 2+ -free/NMDG + -containing medium. S2 Fig shows that in the absence of Ca 2+ , a 40-s application of ATP (100 μM) at −60 mV in bi-ionic NMDG + out/Na + in solution (where the reversal potential of ATP-induced current is about -70 mV [27]) evoked only outward Na + current, whereas in the presence of IVM, outward Na + current was followed by inward NMDG + current. These experiments do not argue against findings with P2X2R presented recently [13] but provide evidence for the existence of two conductive pore states of P2X4R. These pore states, termed open 1 and open 2 , differ in their selectivity for organic cations (a consequence of altered relative permeability P NMDG /P Na ), and the priming of receptors by IVM is needed to switch from one conductive state to another.
We will demonstrate later that the shift in E rev does not require an increase in unitary conductance associated with the open 2 state(s), but rather depends on the selectivity associated with Na + and NMDG + , as suggested by the Goldman-Hodgkin-Katz equation, where R is the gas constant, T is the absolute temperature, F is Faraday's constant, P Na (P NMDG ) is Na + (NMDG + ) permeability, and [Na + ] out ([Na + ] in ) is Na + concentration outside (inside) the cell, whereas [NMDG + ] out ([NMDG + ] in ) is NMDG + concentration outside (inside) the cell. Because the experimentally observed E rev shift is independent of the increase in unitary conductance, the term open 2 state will be used to refer to both the (small) increase in unitary conductance and the (large) change in ion selectivity of the P2X4R pore.

Desensitization masks the increase in unitary conductance in P2X4R
The previous paragraph proposed that P2X4R opens with the open 2 pore state(s) in the presence of IVM, which may have an increase in unitary conductance of as much as 20% [17]. At the same time, the ramp protocol shows a decrease in the slopes of the I-V curves (S1 Fig).
To enforce such an outcome in any potential model of P2X4R  [19]. We are thus led to assume that the probability of finding open receptors on the cell membrane, P(open 1 ), is a strictly decreasing function of time ð _ Pðopen 1 Þ < 0Þ. On the other hand, we expect that the probability of finding a receptor whose pore is in the open 2 state(s), P(open 2 ) to be an increasing function of time ð _ Pðopen 2 Þ > 0Þ. Without specifying a Markov model to describe P2X4R kinetics, we may consider a generic equation for current production in these receptors, capable of distinguishing open 1 and open 2 states based on their conductances and reversal potentials established after washout of intracellular Na + and gain of NMDG + . According to the description above, we can write the equation for current as where g 1 is the maximum conductance of the open 1 state(s), g 2 (> g 1 ) is the maximum conductance of open 2 state(s), and E 1 and E 2 are the reversal potentials associated with the open 1 and open 2 states, respectively. The current equation can be rewritten in a standard form to isolate the total conductance and reversal potential of the cell, as follows where g tot and E tot are the total conductance and reversal potential of the cell, respectively. By equating Eqs (1) and (2), we obtain The requirement for the slope of the I-V curves to decrease during the ramp protocol can be met if the total conductance of the receptor population decreases over time, i.e., _ g tot < 0. Taking the time derivative of g tot and rearranging the terms, we obtain which is strictly negative if we impose the condition It follows that As a first approximation, we can attribute the decrease in the fraction of open states to two processes, desensitization and priming of receptors, related by the equation where δ is the rate of change of open receptors due to desensitization. Furthermore, letting g 2 satisfy where f is the fractional increase in unitary conductance, we can substitute this expression into Eq (3) to obtain Inequality Eq (5) represents a new condition that can be used to produce the decrease in total conductance seen in the ramp protocol. For example, if we consider the experimentally observed value of~0.2 for f in human P2X4R, then the rate of desensitization only needs to be one fifth the rate of the IVM-induced unitary conductance increase in order to mask its effect on the slopes of the I-V curves.
The desensitization rate of naïve receptors is well characterized by the current recordings produced during prolonged application of ATP, which can be used to constrain δ as a fixed parameter. The IVM-induced increase in unitary conductance has not been determined for rat P2X4R, nor its time-course. We determine these in the Markov model (discussed below) by fitting the total current, imposing Inequality Eq (5) to ensure that the total conductance of the cell decreases during the ramp protocol (due to desensitization).

IVM affects receptor desensitization and deactivation kinetics
We next consider the effects of IVM on desensitization and deactivation. During the pulse protocol (Fig 1A), where cells were repeatedly stimulated by 1 μM ATP for 2 s twice per min in the absence (black trace) and presence of 1 μM IVM (colored traces), we observed an initial increase in desensitization rate of the receptor (blue trace in Fig 1A), followed by a gradual decrease in desensitization rate at each subsequent ATP pulse (see the Methods section for quantification procedure). By the fifth pulse (green trace in Fig 1A), the desensitization rate reverted back to a value comparable to that seen in the absence of IVM (black trace in Fig 1A).
To assess if this phenomenon occurs consistently, we evaluated the statistical significance of the transient increase in desensitization rate. To quantify the amount of desensitization seen in the recordings, we used linear fitting to measure the rate of receptor desensitization normalized by the current amplitude of each pulse. As shown in Fig 1B, we did not see a significant change in the rate of desensitization at each ATP pulse in the absence of IVM (filled circles) (n = 7), suggesting that the desensitization proccess of receptors is far from equilibrium. However, in the presence of IVM, the first two ATP pulses following IVM application (indicated by the small arrows) exhibited a significant (p < 0.005 and p < 0.05; n = 7) increase in desensitization rates. The desensitization rate of current recordings in subsequent ATP pulses gradually drifted back to its original value before IVM was applied, further suggesting that the open state, exhibiting an increased desensitization rate, has reached an equilibrium with its corresponding desensitized state. At higher IVM concentrations, however, these transient effects were not observed, but an increase in non-desensitized current amplitude was found [20]. Thus, while the binding of IVM potentiates P2X4R, it also increases both the apparent rate of desensitization at low ATP concentrations and the rate of recovery from desensitization (i.e., it lowers the Gibbs free energy barrier for these transitions).
To assess the deactivation kinetics (i.e., decay of current amplitude following washout of agonist) of P2X4R, we used the same pulse protocol of 1 μM ATP for 2 s twice per minute (Fig 1C). In the absence of IVM (filled circles), receptors underwent fast deactivation with a time constant that remained roughly the same at each pulse, whereas in the presence of 1 μM IVM (open circles following the small arrow) the deactivation time constant progressively increased with incubation time, indicating a decrease in receptor deactivation rates. This effect became even more pronounced at higher IVM concentrations. At IVM concentrations greater than or equal to 10 μM, deactivation following washout of IVM was not always complete (see Fig 2A in [20]), suggesting that complex physiological processes might be initiated at these concentrations. These results are consistent with the idea that IVM increases the sensitivity of the receptor to ATP and decreases the rate of agonist unbinding following its washout from medium [16][17][18].

IVM causes receptor priming during ATP stimulation
Our previous study showed concentration response curves for rat P2X4R stimulated by ATP in the presence and absence of 3 μM IVM (S3 Fig); ATP alone was found to produce a concentration response curve with an EC 50 of 2.3±0.4 μM (blue line), and with 30-s pretreatment with IVM, the EC 50 was 0.5±0.1 μM (green line) [20]. A similar conclusion was reached with human P2X4R [17]. A pretreatment period of 10 s was also considered and, whereas it did produce the same maximal current amplitude, an intermediate EC 50 of 1.6±0.3 μM was measured (maroon line) [20]. This suggests that there are at least two distinct priming effects associated with IVM with separate time scales of action. First, IVM primes receptors by increasing the maximal whole-cell current response. Second, after prolonged exposures to IVM, receptors  [20]). (B) Desensitization rates determined by linear fitting of the desensitization phase of each current pulse. In the absence of IVM (full circles), there was no significant variation in the desensitization rate, but following 1 μM IVM application (arrow), the desensitization rate significantly increased, followed by a gradual decline back to its value prior to stimulation with IVM. Statistical significance was calculated with the one sided Wilcoxon signed rank test. *p<0.05; ***p<0.005. Extended Markov model of P2X4R become further primed by an increased sensitivity to ATP (previously called sensitization). The model in [20] was only partially able to account for these behaviors and specifically was not able to account for their dependence on the duration of pre-treatment because there were no kinetics associated with IVM binding in the absence of ATP.

IVM increases open probabilities
The concentration response curves of P2X4R (S3 Fig) reveal that not only do 10-and 30-s pretreatments with IVM increase sensitivity to ATP (maroon and green lines, respectively), but they also increase the maximum current amplitude evoked by ATP [17,20]. The two hypotheses that can explain this behavior are: (i) the unitary conductance of individual channels increases; or (ii) the number of open receptors is rising (i.e., the maximal open probability increases). Although there is evidence that the former hypothesis holds [17], this does not preclude a change as well in the maximal open probability with IVM application [17,20,23]. In fact, it was reported that the maximal open probability in the absence of IVM is~0.2 compared to~0.8 in the presence of IVM [17,23]. This phenomenon was previously explained by the Markov model in [20], which assumes that IVM modifies the connectivity between open and desensitized states, but that model required a large increase in unitary conductance. This assumption on the conductance is inconsistent with a study of human P2X4R [17], which showed that IVM produces a roughly 5-fold increase in maximal current amplitude while only inducing a 20% increase in unitary conductance. Those authors posited that, in the absence of IVM, desensitization plays a large role in reducing the current amplitude, whereas when IVM is applied, desensitization is greatly reduced, enhancing the observed current. While this is a plausible explanation, we are not aware of any receptors that function in this manner. Moreover, no quantitative analysis was made to assess to what extent such a mechanism produces the observed effect.
In order to test this hypothesis, we constructed a simple and generic mathematical scheme (hereafter referred to as a gating scheme) of a desensitizing ligand-activated receptor (Fig 2A). It consists of two rows: a naïve row comprised of two closed states (C 1 , C 2 ) and two conducting states (Q 1 , Q 2 ), and a desensitized row comprised of four nonconducting desensitized states (D 1 , D 2 , D 3 , D 4 ). As was done in [20], we assume that channels open from states with two or three bound ATP molecules. This is in accordance with the finding that a single-bound receptor state does not lead to activation of P2X7 channels [28]. This is also consistent with previous models of P2XRs and the notion that a single kinetic model underlies the functioning of all receptor subtypes [19,20,[24][25][26]29]. Forward (backward) transitions between two states along each row represent a single ATP binding (unbinding) with rates k 2 , k 4 , k 6 (k 1 , k 3 , k 5 ), respectively, whereas upward (downward) transitions represent desensitization (recovery) with a rate k d (k r ). Concentration response curves were generated for this gating scheme, each with a progressively increasing rate of desensitization k d (see Figs 2A and S5A). It was found that although reduced desensitization rates are capable of increasing the current amplitude at a given agonist (such as ATP for P2X4R) concentration, the mechanism proposed in [17] is unable to significantly increase the maximal current amplitude evoked by the agonist. Rather, it shifts the EC 50 of the concentration response curves leftward as well as increases the Hill coefficient in such a way that the saturating phase of the concentration response curves are shifted by many orders of magnitude. A leftward shift in EC 50 and modulation of the Hill coefficient by IVM have been observed experimentally [17,20]. This was, however, consistently associated with an increase in the maximal current amplitude, which the desensitization mechanism cannot produce at saturating agonist concentrations (see I max in the legend of S5A Fig). A mathematical model introduced by Silberberg et al. also used an IVM-dependent transition rather than modulation of desensitization by IVM in order to produce the increase in maximal current [23]. Therefore, the mechanism suggested in [17] seems unable, at least on its own, to explain the effects of IVM on the concentration-response relationship for the peak current of P2X4R.
After having ruled out decreased desensitization as a cause for the increased maximal current amplitude in the presence of a modulator, we tested an alternative hypothesis, that the closed states exist in equilibrium with a deeply inactivated state (C 0 ) for which the agonist is not effective (Fig 2B). This mechanism has previously been used in Markov models of sodium channels [30]. Transitions linking the two states must be slow, but the equilibrium mixture of the closed-inactivated subsystem (C 0 $ C 1 ) establishes an upper bound on the maximal open probability in the absence of IVM, given by even at the highest agonist concentration. In other words, before agonist application, only some fraction of receptors are in C 1 and are susceptible to agonist-induced activation but in the presence of IVM, more can be recruited (from C 0 ) into C 1 . Evidence of such recruitment was first seen during prolonged application of ATP in the absence and presence of IVM [20] and was obtained also by application of IVM to fully desensitized receptors (S4 Fig). To see how effective this mechanism is in producing the observed effects in [17], we tested the gating scheme of Fig  Extended Markov model of P2X4R C 0 and C 1 ), and plotting the concentration response curves (S5B Fig). Increasing H 1 decreased inactivation, which increased the fraction of receptors in C 1 and thus P Max O . Whereas reducing the occupancy of the deeply inactivated state is highly effective at increasing the maximal current amplitude, it does not significantly shift the concentration response curves or alter the Hill coefficients. Therefore, in order to match the experimental findings that IVM pretreatment of P2X4R not only increased maximal current but also shifted the EC 50 leftward (S3 Fig) and increased the Hill coefficients, both reduced desensitization and rescue from a deeply inactivated state seem to be required. A model incorporating both features is described in the next section.

One-layer model
Based on the above considerations, we designed a one-layer Markov state model that describes the full kinetics of ATP and IVM binding to P2X4R and tested it against experimental data. For a detailed description of the model, see S6 Fig, Table A and Appendix A in S1 Text. Briefly, it is a revised version of the model of Zemkova et. al. [20] that now assumes 3 IVM binding sites, that the binding of IVM acts on P2X4R independently of ATP binding, and that IVM can bind to any ATP-bound state, not just the 3-ATP bound naïve state. Sequential binding of IVM causes three stages of receptor priming, depending on number of IVM molecules bound to receptor: primed-1, primed-2, and primed-3. Primed-1 receptors respond to ATP application with increased current amplitude, reflecting increased open probability. Primed-2 receptors exhibit modestly increased unitary conductance for Na + and significantly increased unitary conductance for NMDG + , whereas primed-3 receptors show increased ATP binding affinity. The model also incorporates rescue from the deeply inactivated state by IVM, and therefore has a maximal open probability given by Eq (6) in the absence of IVM.
Although our analysis of this model (and several variations of it) revealed that it possesses many of the necessary ingredients to capture the gating properties of P2X4R and several aspects of its current recordings (S7 and S8 Figs), it includes the implausible assumption that receptors in the primed states must lose all bound IVM molecules in order to desensitize. This assumption led to two major issues in the performance of the model: (i) it did not capture accurately the short timescales of activation and desensitization robustly; and (ii) it produced discrepancies in current amplitudes when compared to experimental data during the pulse protocol. That motivated us to design a more accurate model of P2X4R kinetics.

Two-layer model
Receptor stages and gating schemes. The observed transient increase in desensitization of P2X4R immediately after IVM application, followed by a gradual decrease in this rate (Fig  1), suggests that IVM bound states of the receptors are also directly associated with desensitization. This type of behaviour cannot be reproduced well by Markov models lacking this feature (such as the model in [20] or the one-layer model). Because it was more reasonable to suggest that IVM modifies receptor desensitization rather than preventing it from happening altogether, we investigated the effect of allowing reversible transitions between IVM-bound states and their corresponding desensitized states. This was done by including, in addition to the IVM-induced modified ATP kinetics (which produce the slowed deactivation rates and increased sensitivity to ATP), another IVM bound row with modified desensitization kinetics that can readily recover from desensitized states. According to both the simple gating schemes of Fig 2 (taken from [20]) and the one-layer model (S6 Fig), the submodel formed by removing the IVM-bound states allows receptors to return from the desensitized row with a rate that is independent of IVM-related processes. We hypothesized that the gating scheme in Fig 2A is the underlying system whose rate parameters are modified with each IVM binding to produce the varied kinetics we observe at different IVM and ATP concentrations.
To test this, we chose two recordings of the pulse protocol performed at 1 μM and 10 μM IVM (see the experimental recordings in S9 Fig). Choosing the first pulse (before IVM application) and the last pulse (before IVM washout) from each recording, respectively, we were able to fit the transition rate parameters of the gating scheme in Fig 2A to these pulses very accurately (see S9A and S9C Fig). In particular, this model captured the fast activation, desensitization, and insensitivity to ATP removal at 10 μM IVM with one set of parameters for each case. However, S9B Fig shows that although the gating scheme reproduced the first pulse extremely well, by the last pulse, both the deactivation kinetics and current amplitude were off. Similarly, S9D Fig shows good agreement between the last pulse and simulation while the first few pulses of the recording were off. These results suggest that a mixture of these gating schemes must be linked by IVM dependent transition rates, and that (i) at the lowest IVM concentrations, the mixture will be predominantly composed of naïve receptors and receptors which have only undergone a single modification by IVM, and (ii) at the highest IVM concentrations, the mixture will eventually saturate at these states that are most modified by IVM. In other words, we expect that building a model in which gating schemes of the type depicted in Fig 2A are allowed to mix, will better capture the complex experimental behaviour.
Description of the two-layer model. Motivated by the fact that the same underlying gating scheme could be used to reproduce single pulses of the pulse protocol with great fidelity at all IVM concentrations, we set out to develop a model which could capture not only the short timescale behaviour of receptors but also the progression in kinetics and amplitude which the single gating scheme failed to capture. When cooperativity between IVM and ATP was included in simple models (such as the one-layer model), there was a clear progression in kinetics and amplitude. With this in mind, we have investigated a new model that included a mixture of gating schemes shown in Fig 2A and 2B. Because we allow for all closed and open states to desensitize, we termed this model the two-layer model. Fig 3A describes the two-layer model in detail. P2X4Rs are assumed to be found in four stages during their activation cycle: deeply inactivated, functional, desensitized, and internalized. The upper layer consists of closed (C 1 −C 8 ) and conducting (Q 1 −Q 8 ) states, which we termed functional states. Fig 3B shows how the functional states are determined by binding of ATP and IVM, producing four subgroups of functional states: naïve, primed-1, primed-2, and primed-3. Naïve and primed-1 channels, which differ in their maximum probability of opening, respond to binding of 2 or 3 ATP molecules with open 1 states. Primed-2 and primed-3 channels respond to ATP with open 2 (NMDG + -conducting) states, but primed-3 receptors also exhibit increased sensitivity to ATP. The model was constructed by linking 4 distinct gating schemes together (each consisting of 4 functional states and 4 desensitized states) through fourth-order Hill functions that depend on IVM concentrations (see Appendix B and Table C in S1 Text). Their lower layer is composed of nonconducting desensitized states. This effectively means that desensitization is a process which does not modify agonist binding but rather affects the receptor-channel by making it nonconducting despite agonist binding. Furthermore, we allowed for internalization from all desensitized states which have three bound ATP molecules. The rate of internalization from the naïve states (H 3 ) was assumed to be different from that of IVM bound states (H 5 ) to reflect that IVM alters endocytosis of P2X4R receptors [31].
The results of Fig 1 imply that IVM not only modifies ATP-binding and activation kinetics, but also shifts the receptor population into states where the transition to desensitized states is less favorable. To obtain such behaviour in the model, we must ensure that the ratio of desensitization-to-recovery rates associated with the naïve states is much larger than those of the primed states, i.e., the following condition must be satisfied: Fitting of the two-layer model to the pulse/prolonged protocols. In order to assess the capacity of the model to reproduce P2X4R gating, we used MCMC techniques (see Methods Section) to determine the parameter values that reproduce the kinetics observed experimentally during the pulse and prolonged protocols. The results shown in Figs 4 and 5 reveal that the model is able to capture every aspect of P2X4R gating, including naïve receptor activation, desensitization, and deactivation ( Fig 4A, 4E and 4F) in the absence of IVM, and the increase in current amplitude (Fig 4B-4E) and deactivation time constant (Fig 4F) upon IVM application during the pulse protocol (compare experimental recordings in blue curves/white bars to simulations in red curves/black bars). The transient increase in desensitization rate after 1 and 10 μM IVM application followed by a progressive decrease in desensitization rate at each subsequent pulse (Fig 4G and 4H), along with the insensitivity of P2X4R to ATP removal during washout (Fig 4H) in the pulse protocol were also captured by the model. The ATP-dependent primed with 1, 2 or 3 molecules of IVM. States Q 1 ,Q 2 ,Q 3 and Q 4 are open with maximum conductance g 1 and reversal potential E 1 (open 1 ) whereas Q 5 ,Q 6 ,Q 7 and Q 8 are open with maximum conductance g 2 > g 1 and reversal potential E 2 < E 1 (open 2 ). Before stimulation with ATP and IVM, receptors reside in C 0 or C 1 . Parameter values are listed in Table C   concentration response curves for current amplitude generated by the model during the pulse protocol in the absence (black) and presence (red) of IVM applied 30 s before ATP stimulation (Fig 5A) also exhibited a leftward shift (from~1.78 μM to~0.587 μM) in the EC 50 in a manner almost identical to that seen experimentally (dotted lines). These results suggest that the model can reproduce the kinetics of the receptors occurring at short time scales, that having desensitization decoupled from IVM unbinding allowed the model to produce the transient increase in desensitization during low IVM stimulation (i.e., in primed-1 cells) and that (pre)stimulation with moderate to high IVM concentrations rescues receptors from desensitization (through an increase in the rate of return from desensitization) within this time scale.
Comparing the observed data during the ramp protocol to the outcomes of the model shows that the model can reproduce the decrease in the slope of the I-V curves as demonstrated by the first and last I-V curves generated by applying a voltage ramp from − 80 mV to + 80 mV during 10 s of 100 μM ATP stimulation in the presence of 3 μM IVM applied 20 s prior to ATP stimulation (Fig 5B). Here we focused on the first and last I-V curves because we are not tracking ionic fluxes within the model (see Appendix C in S1 Text). This phenomenon was robustly produced by ensuring that Inequality Eq (5) was satisfied by model simulations. The current equations (Eq (S1)) can also produce a shift in reversal potential from −36.7 mV to − 27.5 mV, similar to that observed experimentally in the presence IVM under bionic (NMDG + out/ Na + in) conditions. Based on the Goldman-Hodgkin-Katz equation, this was achieved by generating a value for E 1 , the reversal potential of the open state (open 1 ), using the MCMC algorithm, that is more negative than that of the primed state E 2 (see Appendix C in S1 Text). It is important to note that (i) this shift in reversal potential can be still generated in the absence of an increase in unitary conductance, when g 1 = g 2 , as shown in S12B Fig, indicat-ing that the shift is independent of an increase in unitary conductance; and (ii) we do not explicitly model the flux of cations through the pore, and therefore the reversal potential E 2 reflects both the altered selectivity of the open 2 pore and the accumulation of NMDG + within the cell (as was demonstrated in [13]).
Finally, comparing the experimental recordings (blue) to model simulations (red), shows that the two-layer model is successful in reproducing the responses obtained in prolonged application of 100 μM of ATP in both the absence (Fig 5C) and presence of 3 μM IVM ( Fig  5D). The former showed partial recovery of peak current response after 4-min ATP washout, whereas the latter showed a significant increase in amplitude and a plateauing in current amplitude towards the end of the desensitization phase, as seen experimentally. These results suggest that a substantial fraction of P2X4R receptors are shielded from desensitization in the presence of IVM during prolonged time scales (by transitioning receptors to the primed gating schemes) and that these four processes exist at equilibrium during the plateauing phase.
Kinetics of the two-layer model. The structure of the two-layer model and the fitted parameter values (Table C and Appendix B in S1 Text) suggest that agonist potency is increased by IVM which rescues receptor population from a pool of inactivated receptors (C 0 ). The gating scheme of the primed-1 subsystem (i.e., with one IVM bound) had an increased desensitization rate which led to the observed transient increase in desensitization after IVM application and was also found to exhibit negative cooperativity (see Methods). The increased sensitivity to ATP in this and other IVM-bound schemes was produced by having their ATP unbinding rates decrease and/or their ATP binding kinetics increase. The model also revealed that the recovery rates of all the IVM-bound gating schemes are significantly larger than that of the naïve gating scheme (i.e., k r,i ) k r,1 , i = 2, 3, 4) while the internalization rate of these schemes is lower than that of naïve receptors (i.e., H 5 < H 3 ), which leads to the long-lasting currents observed during prolonged ATP applications in the presence of IVM (Fig 5D; for additional current recordings with the same protocol, see S10 Fig). Finally, the unitary conductance increase estimated by MCMC was approximately 15%, which is within the experimental bound. Together, these results show that, in this model of P2X4R, unitary conductance increase is present but plays a minor role in increasing the current amplitude in the presence of IVM, whereas the open probability is the key to producing such large currents.

Discussion
Data-based Markov state models that describe the processes of ligand binding/unbinding to ligand-gated receptor are powerful tools to understand orthosteric and allosteric regulation of these channels. P2X4Rs are prototypical examples of such receptors with orthosteric and allosteric binding sites for ATP and IVM, respectively. They are associated with ion channels that are permeable to small cations, including Na + , Ca 2+ , Mg 2+ and K + . The binding of ATP leads to receptor activation and channel opening, while IVM binding increases receptor unitary conductance and sensitivity to ATP. Here we analysed the kinetics of ATP and IVM binding/ unbinding to P2X4R, and determined its gating properties using two detailed Markov models labelled the one-layer and two-layer models.
The one-layer model extended a previously developed, simple Markov model of P2X4R by taking into consideration a deeply inactivated state, nonresponsive to ATP but responsive to IVM (existing in equilibrium with the naïve ATP-unbound state), along with three additional gating schemes (per each ATP binding) representing the three IVM binding sites. The model also assumed that the IVM and ATP binding are independent of one another and that sequential binding of IVM can occur at any ATP-bound or unbound states. Our analysis revealed that the deeply inactivated state was essential for capturing the increase in the maximum response (I max ) in the ATP-dependent concentration-response curves (in the presence of IVM), with only a small increase in conductance between the open 1 and open 2 states. Although the model was able to capture many of the essential features of P2X4R recordings (S7 and S8 Figs), it assumed that IVM bound states can only desensitize by first becoming completely free of IVM. That made the model unable to robustly capture the short timescales of activation and desensitization, and it also produced current amplitudes incompatible with experimental data during the pulse protocol.
By allowing the IVM-bound states to desensitize, we were able to show that a simple gating scheme is able to capture the profiles of the first pulse (in the absence of IVM) or the last pulse (in the presence of IVM) of the pulse protocol very accurately when the scheme is fitted individually to each pulse, but not both simultaneously. However, when comparing the entire pulse protocol recording to the outcome of the scheme for each case, there was a gradual increase in discrepancy between them, suggesting that a mixture of gating schemes must coexist to be able to capture all aspects of P2X4R kinetics. That led us to propose the two-layer model, which assumes that IVM-bound states can desensitize.
The two-layer model was successful in capturing every aspect of P2X4R kinetics very accurately, including the short and long time scales of activation and desensitization, particularly the changes in the desensitization rate observed during the pulse protocol of Fig 1, as well as the current amplitudes. The observed shift in the EC 50 along with the increase in the maximum current, during pre-stimulation with IVM, were also reproduced by the model (through the presence of the deeply inactivated state). Moreover, these gating schemes can be used to understand why ATP binding mutants with low amplitude of response tend to have significantly larger fold-increases in maximal current in the presence of IVM [21]. If we view such mutants as disproportionately populating the deeply inactivated state, where they cannot bind ATP, then their rescue by IVM from this state will produce a much larger fold-increase in maximal current. The existence of this deeply inactivated state was probed experimentally, by applying IVM to a cell whose receptors were almost completely desensitized from prolonged applications of 100 μM ATP. Upon IVM application, we observed an increase in the maximal current amplitude to about half of the initial maximal current (see S4 Fig), suggesting that IVM rescued receptors from a (deeply inactivated) pool corresponding to about one third of all receptors.
The two-layer model consisted of 4 gating schemes linked together through ATP/IVM binding/unbinding. Two of these gating schemes (the primed-2 and -3) contained conducting states exhibiting a 15% increase in unitary conductance compared to that of the open states in the naïve and primed-1 states. This increase is within the 20% limit seen experimentally [17], and is not required to produce IVM's increase in maximal current amplitude within our model. Instead, the effect of IVM on maximal current amplitude is produced mainly by an increase in open probability. According to Eq (6), the two layer model predicts a maximal open probability of approximately 0.53 in the absence of IVM (with 47% of receptors in a deeply inactivated state), while it can easily reach values greater than 0.9 in the presence of IVM. It was suggested in [13] that the ionic conditions in the medium and the pipete are responsible for producing electrochemical effects which were long presumed to be evidence for pore dilation, particularly in P2X2R. While temporal changes in ionic gradients play a significant role in producing a shift in E rev associated with the I-V curves during the ramp protocol, our results suggest that the transition to open 2 (which is permeable to NMDG + ) is an intrinsic property of the pore in P2X4R, is independent of the increase in unitary conductance (S12B Fig) and is induced by priming with IVM. The two-layer model assumes that the increase in unitary conductance associated with this transition is masked by desensitization. This results in the shift in E rev being accompanied by a decrease in the slope of the I-V curves (due to desensitization). The observed decline in the slope of the I-V curves was a consequence of Inequalities Eqs (3) and (5). The previously developed model in [20] was capable of fulfilling these conditions and producing the decrease in the slope of the I-V curves, but it required a large increase (>150%) in unitary conductance to achieve it while simultaneously producing the increase in current amplitude induced by IVM. The inclusion of the primed-1 row with conductance g 1 (and reversal potential E 1 ) was an essential element for reproducing the shift in the I-V curves. Without this intermediate step, the model required a very positive reversal potential for the open 2 state (E 2 ), which does not reflect its loss of selectivity. Both the onelayer and two-layer models proposed here keep the increase in unitary conductance within 20% and produce the current growth with IVM pretreatment through IVM-induced transitions from the inactivated state C 0 (by increasing P Max O , as given by Eq (6), rather than increasing the maximal conductance). The two-layer model, however, is more plausible because it does not assume that desensitization necessitates the unbinding of IVM. Moreover, according to this model, IVM is able to transition receptors to the primed-3 states in the absence of ATP, allowing pretreatment with IVM to produce sensitization independently of ATP. It is important to note that effectively removing the 15% increase in unitary conductance associated with the primed states only slightly altered model simulations and did not abolish any experimental phenomena in symmetric ionic conditions (S11 and S12 Figs). Thus, the results of the twolayer model are independent of an increase in unitary conductance, but require a change in selectivity in order to capture the shift in the I-V curves of the ramp protocol.
Recently, the pore dilation hypothesis has become increasingly disputed. Molecular dynamics simulations indicate that NDMG + is capable of permeating the open state of P2X4R pore, provided it is maintained in an open state long enough for the slow permeation event of NDMG + to take place [15,32]. One of the primary effects of IVM application on the single channel kinetics of P2X4R is to shift the distribution of open times from the sub-millisecond timecale to tens of millisecond [17]. We hypothesize that the drastic change in P2X4R's permeability for NMDG + upon IVM application results from the priming of receptors in such a way that their pores remain in the open 2 state for long enough for the slow permeation of NMDG + . This increase in the permeability for NMDG + (via an IVM-dependent transition to the open 2 state) allows for its influx into the cell. Together with the efflux of Na + , these fluxes produce a more positive E rev as determined by the Goldman-Hodgkin-Katz equation.
While the two-layer model may seem to be a large departure from both the previously developed model in [20] and the one-layer model, it should be noted that in the absence of IVM, the remaining blocks of the models (or submodels) are identical, and that the P2X2R Markov model developed in [19] has a similar structure; it included a corresponding desensitized state for each of its closed and open states, although the desensitization pathway for primed (sensitized) states was calcium dependent. The increase in the number of states and number of kinetics parameters in the two-layer model was necessary to capture all the observed features of P2X4R, which previous models, including the one-layer model presented here, failed to do. A step by step validation of such an increase in complexity was provided through the use of coupled gating schemes, and the design of an extensive MCMC fitting algorithm that combined parallel tempering approaches with the t-walk method to estimate the kinetic parameters of the model efficiently.
The two major allosteric effects of IVM's on P2X4Rs, observable from whole cell currents as an increase in maximal current amplitude or the deactivation time constant, exhibit distinct concentration dependencies [17,20]. This suggests that they are likely caused by two independent processes. This existence of two distinct allosteric effects with differing concentration dependencies have also been reported for other P2XRs [11,19], although for these receptors, ATP alone was sufficient to induce such effects. The models presented here reduce the concentration dependence of IVM's allosteric effects on P2X4R to a single sequential binding process, and thus represent a major simplification of a more realistic model where all effects are assumed to arise from independent binding events. Despite this simplification, the model is quite capable of capturing all aspects of allosteric modification by IVM.
An important item for future work is cooperativity in the ATP and IVM binding, which has been investigated in the one-layer model but not yet in the two-layer model. By assuming correlations between the binding/unbinding parameters of ATP and IVM in the one-layer model, we were able to reduce the number of estimated parameters in that model significantly and found that there was negative cooperativity in the ATP binding in the naïve and primed-1 rows. Investigating if such cooperativity exists in the two-layer model is also warranted. This can be done by imposing correlations between the kinetic parameters of the two-layer model, which will again reduce the number of estimated parameters, and testing for cooperativity in ATP binding, IVM binding and between ATP and IVM binding. These variations of the model can then be compared to each other using Bayesian approaches to determine which is most likely.
In conclusion, here we present two novel models, one of which (the two-layer model) effectively mimics all experimental observations. In this model, receptors go through four stages of activation cycle during ATP and IVM binding: transitioning from functional to desensitized, from desensitized to internalized, from internalized to deeply inactivated and from deeply inactivated to functional. Functional and desensitized stages each exist as 16 distinct states, determined by the progressive saturation of three ATP and three IVM binding sites, whereas internalized and deeply inactivated receptors are single states. Binding of IVM influences ATP-induced gating properties of receptors, i.e. the rates of activation, desensitization and deactivation, open probability of channels, and the sensitivity of receptors to ATP. The channel pore state, open 1 is predominantly permeable to small cations and open 2 is permeable to large organic cations.

Cell cultures and transfection
Experiments were performed on human embryonic kidney 293 cells (HEK293; American Type Culture Collection), which were grown in Dulbecco's modified Eagle's medium supplemented with 10% fetal bovine serum, 50 U/mL penicillin, and 50 μg/mL streptomycin in a humidified 5% CO 2 atmosphere at 37˚C. Cells were cultured in 75-cm 2 plastic culture flasks (NUNC, Rochester, NY, USA) for 36-72 h until they reached 80-95% confluence. Before the day of transfection,~150 000 cells were plated on 35 mm culture dishes (Sarstedt, Newton, NC, USA) and incubated at 37˚C for at least 24 h. For each culture dish of HEK293 cells, transfection of wild-type P2X4R was conducted using 2 μg of DNA with 2 μl of jetPRIME reagent in 2 ml of Dulbecco modified Eagle's medium, according to the manufacturer's instructions (PolyPlus-transfection, Illkirch, France). After 24-48 h of incubation, the transfected cells were mechanically dispersed and re-cultured on 35 mm dishes of Corning 3294 CellBIND Surface for 1-4 hours before recording. Transfected cells were identified by the fluorescence signal of enhanced green fluorescent protein using an inverted research microscope with fluorescence illuminators (Model IX71; Olympus, Melville, NY).

Patch-clamp recordings
Currents were recorded in a whole-cell configuration from cells clamped to −60 mV using an Axopatch 200B patch clamp amplifier (Axon Instruments, Union City, CA, USA). All currents were captured and stored using Digidata 1550A and pClamp10 software package. Patch electrodes were pulled from borosilicate glass tube with a 1.65 mm outer diameter (type GB150F-8P; Science Products GmbG, Hofheim, Germany) using a Flaming Brown horizontal puller (P-97; Sutter Instruments, Novato, CA). The tip of the pipette was heat-polished to a final tip resistance of 3-5 MOhm. During the experiments, the dishes with cell cultures were perfused with an extracellular solution containing: 142 mM NaCl, 3 mM KCl, 2 mM CaCl 2 , 1 mM MgCl 2 , 10 mM HEPES and 10 mM D-glucose, adjusted to pH 7.3 with 10 M NaOH. The osmolarity of solution was 290 mOsm as determined by a vapor pressure osmometer (Model VAPRO 5520; Wescor, Logan, UT, USA). Experiments were done on single cells with an average capacitance of about 10 pF held at membrane potential of -60 mV. Patch electrodes used for whole-cell recording were filled with an intracellular solution containing: 145 mM NaCl, 10 mM EGTA and 10 mM HEPES; the pH was adjusted with 10 M NaOH to 7.

Mathematical models
The one-layer Markov state model developed here (S6 Fig), is a revised version of a previously developed Markov model describing P2X4R orthosteric activation by ATP and allosteric modulation by IVM [20], whereas the two-layer model (Fig 3) is more complex in nature, taking into consideration all processes involved in ATP and IVM binding. Both models assume the presence of three ATP and three IVM binding sites and were tested against current recordings to compare their performance in capturing the physiological properties of P2X4R. The following symbols were used to describe the various states of the model: C for closed, Q for conducting (open 1 and open 2 ), D for desensitized and Z for internalized states, each representing the fraction of receptors in a given state. The transition rates between the various states are in Tables A and C in S1 Text. Detailed balance was not explicitly incorporated into these models because it was assumed that P2X4R never reach absolute equilibrium during ATP and IVM stimulation. Stimulation with ATP during the pulse protocol (repetitive stimulation with 1 μM ATP for 2 s twice per minute in the absence and presence of various concentrations of IVM applied to the bath medium after two ATP pulses) and the prolonged protocol (stimulation with 100 μM ATP for extended periods, longer than 1 min, in the absence and presence of 3 μM IVM) were modeled as a square wave and a rectangular function for the duration of ATP application, respectively; stimulation with IVM was modeled as a rectangular function for the duration of IVM application. The ramp protocol was modeled as a sawtooth-like sequence of upstrokes with slope 320 mV/s (rising from −80 mV +80 mV over 500 ms). Detailed descriptions of the two Markov models along with their differential equations are provided in Appendices A and B in S1 Text.

Data characterization
Concentration-response data points were fit to a hill function where y is the amplitude of the current evoked at a particular ATP concentration x, I max is the maximum current observed at 100 μM ATP, EC 50 is the ATP concentration producing 50% of the maximum current, and n is the Hill coefficient. Deactivation kinetics of the current decay after agonist washout were fitted to a single exponential or to a sum of two exponentials where A 1 and A 2 are the amplitudes of decay for the first and second exponentials, τ 1 and τ 2 are their decay time constants, and C is the baseline current. In the case where the sum of exponentials fits the data better than a single exponential, we report the weighted time constant In either case, we labeled the derived time constant of deactivation as τ off . Statistical significance ( ÃÃ p<0.01 and Ã p<0.05) was assessed using the Wilcoxon signed rank test. MATLAB (MathWorks, Natick, MA) was used to solve the differential equations of the models numerically, fit the models to the data and apply statistical tests.
In order to quantify the rate of desensitization from current traces of the pulse protocol (see Figs 1 and 4) which do not show complete desensitization, we note that, at low concentrations of ATP (1 μM), desensitization is a mono-exponential process and thus employ a mono-exponential model where A is the magnitude of the current at the onset of desensitization and τ is the time constant of desensitization. By first normalizing the current and then evaluating its derivative, with respect to time, at the onset of desensitization (t = 0), we obtain Thus the first derivative at time t = 0 of each desensitizing current, _ Ið0Þ, normalized by its maximum current (as has been plotted in Fig 1B), yields information about the time constant of desensitization (i.e., a small normalized initial desensitization rate corresponds to a large time constant of desensitization). Due to the simultaneous activation and desensitization of multiple receptors, _ I ð0Þ was estimated using a linear fit for the small window of time (1-1.5 s) after cells achieved their maximal current and before agonist was removed. This window is relatively small compared to the 6 s time constant of desensitization for P2X4R, and therefore this approximation method provides a relevant estimate of _ I ð0Þ=A, which serves as a proxy for desensitization time constant τ.

Parameter estimation (Adaptive Parallel Tempering & t-Walk)
Parameter estimation was performed using MCMC techniques. Model simulations were generated using ode solvers in MATLAB and then fit to experimental recordings. Generally, MCMC produces Markov chains Λ = {x 1, where N is the number of parameters (p i ) and m = 1, 2, Á Á Á, M is the m th iterate of the Markov chain. The iterates represent samples from the posterior distribution π(x) determined using Bayes' theorem as follows pðxÞ / LðxÞPðxÞ; where L(x) = P(data | x), the likelihood function, is the probability of observing the data given the parameter values of x, and P(x) is the prior distribution of x, which reflects any prior knowledge about the parameter values independent of observed data. Proportionality, indicated by /, is sufficient-therefore there is no need to normalize the posterior.
In order to increase mixing of modes in parameter space, we used the parallel tempering algorithm which produces Markov chains in the product space where each chain x ðlÞ m , l = 1, 2, Á Á Á, L, was sampled from a tempered distribution p b ðlÞ and β (l) is the inverse temperature of each chain. Parameter sets were stochastically swapped between chains according to the swap kernel of Miasojedow et al., and their strategy of adaptively updating the inverse temperature of each chain [33] was adopted. Because Metropolis-Hasting move-kernels can be difficult to tune for continuous-time Markov models of ion channels [34], we used the adaptive move kernel of the t-walk sampler [35] instead. Since the t-walk samples from the product distribution π(x)π(x 0 ), the composite MCMC method samples from the product distribution Given that we have a set of discretely sampled whole-cell current recordings, we initially adopted the likelihood function from Gregory [36], defined by where the index i refers to the i th discretely sampled data point, K is the number of data points in the experimental recording, I μ,i and σ i are the i th samples, respectively, of the mean current and current standard deviation estimated from the data set, and I x,i is the i th sample of the current produced by the model given the set of parameters x. To circumvent (i) sampling inefficiency from the posterior distribution, which is exacerbated by the use of high data sampling rates, and (ii) very poor fitting of rapid transient behaviour, due to the value of the likelihood being dominated by slower portions of the signal with more data points, we opted to fit (using the least-squares method) both experimental and simulation data to appropriately chosen functions and to compare the fit parameters of the experimental and simulated data. For example, we have used exponential functions to measure deactivation kinetics (as described above). where τ off , μ and σ τ are the mean deactivation time constant and its variance, respectively, A off,μ and σ τ are the mean deactivation amplitude and its variance, respectively, and τ off , x and A off , x are the deactivation time constant and amplitude produced by the model corresponding to the parameter values x.Using this description-based approach (rather than the distance from data), we were able to simultaneously fit numerous aspects of P2X4R activation kinetics and their allosteric modulation. This was done by comparing experimental data and model predictions of (i) time dependence of the activation time and normalized rate of desensitzation in the absence and presence of 1 μm IVM (Fig 1A and 1B) (ii) maximal current, deactivation time constant, and desensitization at 1 μM ATP with increasing IVM concentrations (Fig 4E and  4F) (iii). Insensitivity to ATP removal at 10 μM IVM (Fig 4H) (iv) decay of current amplitude deactivation time constant following IVM washout (Fig 4A-4D) (v) activation, desensitization, and recovery after washout of 100 μM ATP in the absence and presence of IVM (Fig 5C and  5D) (vi) EC 50 and Hill coefficient (n) of the ATP concentration-response curves for peak current ( Fig 5A).
The degree of cooperativity in ATP binding was determined from the Markov chain x m of 5000 samples associated with each ATP binding and unbinding rate k i , i = 1, 2, Á Á Á, 24, that was generated from data fitting, followed by calculating the chains of ATP binding affinities along each of the non-desensitized rows of the one-layer model (including the four lower rows) and the two-layer model (including the four rows in the upper layer). The posterior distributions associated with these affinities were used to compare the values of the most frequently sampled points along each row (n = 1, 2, 3). In the presence of a specific cooperativity between ATP bindings, correlations between the different binding affinities were detected and reported. Hill functions were fit to ATP-dependent concentration-response curves of peak current amplitude generated using the gating schemes in Fig 2A and 2B, respectively. Enhanced desensitization and inactivation were realized by altering the magnitude of k d in A and H 1 in B, respectively, by a factor indicated in the legends. The fitted parameters of each Hill function (including I max , EC 50 and n) for both cases are also presented in the legends. Enhanced desensitization increases EC 50 and n, whereas enhanced inactivation decreases I max .

Supporting information
Note that the range of agonist concentration along the x-axis spans 6 orders of magnitude but all responses in B reach saturation within 2 orders of magnitude, as observed experimentally (compare to S3 Fig There are three ATP and three IVM binding sites; for details see Fig 3B. The states C i , D j and Q i , i = 1, 2, Á Á Á, 8 and j = 1, 2, 3, 4, correspond to the fraction of receptors in the closed, desensitized and conducting states, respectively, and the states C 0 and Z correspond to the fraction of receptors in the deeply inactivated and internalized states, respectively. The states Q 1 , Q 2 , Q 3 and Q 4 are open with maximum conductance g 1 and reversal potential E 1 , whereas the primed states Q 5 , Q 6 , Q 7 and Q 8 are open with maximum conductance g 2 > g 1 and reversal potential E 2 > E 1 . Transitions between states along rows represent binding (rightward) and unbinding (leftward) of one ATP molecule, and along columns (below the naïve row) represent the binding (downward) and unbinding (upward) of one IVM molecule. Before stimulation with ATP and IVM, receptors reside in either C 0 or C 1 . Primed receptors must transition through the naïve states to desensitize, which is unlikely. The two-layer model (Fig 3A) removes this constraint. Parameter values of the one-layer model are listed in Table A