Kinetics of Extracellular ATP from Goldﬁsh Hepatocytes: A Lesson from Mathematical Modeling

In goldﬁsh hepatocytes, hypotonic exposure leads to cell swelling, followed by a compensatory shrinkage termed RVD. It has been previously shown that ATP is accumulated in the extracellular medium of swollen cells in a non-linear fashion, and that extracellular ATP (ATPe) is an essential intermediate to trigger RVD. Thus, to understand how RVD proceeds in goldﬁsh hepatocytes, we developed two mathematical models accounting for the experimental ATPe kinetics reported recently by Pafundo et al. in Am. J. Physiol. 294, R220–R233, 2008. Four different equations for ATPe ﬂuxes were built to account for the release of ATP by lytic ( J L ) and nonlytic mechanisms ( J NL ), ATPe diffusion ( J D ), and ATPe consumption by ectonucleotidases ( J V ). Particular focus was given to J NL , deﬁned as the


Introduction
Maintenance of cell volume in the face of fluctuating intra and extracellular osmolarity is essential for cell function.It is especially significant for the animal cell because they cannot withstand osmotic gradients across their membranes.Thus, it is not surprising to find that cells have developed mechanisms that sense and oppose volume changes.Since in most animal cells net water movements across the cell membrane are driven almost exclusively by the osmotic gradient, a reduction of extracellular osmolarity will lead to fast cell swelling.However, during continuous hyposmotic stress, cells undergo a slower secondary shrinkage known as regulatory volume decrease (RVD), thereby preventing cell disruption, and thus maintaining cell viability.Although the signaling pathways activating this regulatory response have not yet been elucidated fully, RVD is mediated to a large extent by KCl loss through electroneutral ion transport pathways or by the separate activation of K + and anion channels (Haussinger, 1996;Jakab et al., 2002).
Recently, there has been a growing body of evidence showing that extracellular nucleotides, mainly ATP, play a significant role as extracellular signals in cell volume regulation, where they act by binding to specific cell surface P receptors (purinic and pyrimidinic receptors, Lazarowski et al., 2003).At the same time, it was shown that most cells are capable of releasing ATP to the extracellular space, particularly when cells are challenged by hypotonicity.Accordingly, Wang et al. (1996) presented the first experimental evidence to support a model for autocrine and paracrine modulation of RVD by extracellular ATP.Using rat hepatoma cells exposed to hypotonicity, cells first swell and then release ATP in response to increases in cell volume.The resulting endogenous extracellular ATP could bind to subtype 2 P receptors to activate Cl − permeability (the effector mechanism) and RVD.The relevance of such findings has been confirmed using primary cultured hepatocytes from humans (Wang et al., 1996) and fish (Pafundo et al., 2004).
A recent study of our group using goldfish hepatocytes exposed to hypotonic medium provided new insights into this nucleotide dependent mechanism of volume regulation (Pafundo et al., 2008).First, we showed that exogenous ATP triggered RVD in a dose dependent manner.Second, swollen 1 goldfish hepatocytes were able to release significant amounts of ATP by both lytic and nonlytic mechanisms.Once in the extracellular medium, extracellular ATP (ATPe) could decrease by diffusion and/or be degraded by ectonucleotidases present at the surface.In addition, we showed that in these cells a positive feedback mechanism exists, whereby extracellular ATP can activate the release of ATP by a yet unidentified mechanism.Thus, in swollen goldfish hepatocytes, the concentration of ATPe located at the cell surface can be affected by several mechanisms interacting simultaneously in a complex way.This is why a simple, preliminary mathematical model was developed accounting for the time course of ATP accumulation in the extracellular space when goldfish hepatocytes are challenged by hypotonicity (Pafundo et al., 2008).
Given the experimental evidence supporting the importance of ATPe as the main trigger of RVD of goldfish hepatocytes, in the present study, we exhaustively investigated and generalized the previous mathematical model.
This work is structured as follows.The first part of the paper is devoted to a general description of the sources (lytic and nonlytic release of ATP) and depletions (ATP consumption by ecto-ATPase activity at the cell surface and ATPe diffusion) of extracellular concentration of ATP ([ATP] e ) as well as their incorporation in the general mathematical model.This procedure resulted in two models described by linear nonautonomous ordinary differential equations.Subsequently, we propose several functions describing the nonlytic release of ATP and considered the presence or not of ATP diffusion within the extracellular compartment.In order to check if these functions as well as the presence or not of ATP diffusion can produce results with physical meaning, the models are evaluated at infinite time.After that, simulations are performed to obtain the best values of the functions parameters.Finally, the best model is chosen after the running model dependent fits to recently published experimental data of our group (Pafundo et al., 2008).The best model obtained allowed us to quantify the contribution of the different processes that determine the kinetics of [ATP] e and at the same time predict different physiological scenarios that goldfish hepatocytes may face.Moreover, since there is up to date no molecular identification of ATP secretion pathways in liver cells, the model is useful in that it includes several kinetic properties that a putative ATP transporter must fulfill.

Theory
In this section, mathematical models were used to analyze the [ATP] e kinetics of goldfish hepatocytes challenged by hypotonicity.In order to analyze the contribution of ATPe diffusion in the modeling, two main types of models were developed, i.e., a model without (model 1) and a model with (model 2) ATPe diffusion.Both models were built considering several ATP fluxes that either increase or decrease the [ATP] e .

Model 1 (without diffusion)
This is a two-compartment model, involving one intracellular and one extracellular compartment.Diffusion of ATPe within the extracellular medium is considered to take place instantaneously.The [ATP] e is a state variable that can be controlled by the following ATP fluxes (Fig. 1A).
(1) Flux of nonlytic release of ATP (J NL ).That is, a release of ATP from the cells not involving cell membrane rupture.
(2) Flux of lytic release of ATP (J L ).The loss of cell viability leads to the release of ATP to the extracellular medium.
(3) Flux of ATP consumption by ecto-ATPase activity (J V ), to account for ATP hydrolysis at the cell surface.2.1.1.Properties of the ATPe fluxes 2.1.1.1.Nonlytic ATP flux Viable goldfish hepatocytes undergo a nonlytic release of ATP (J NL ) in response to a hypotonic challenge.Since the molecular mechanism mediating this flux has not been identified, we evaluated a set of trial time functions (collectively termed as J R (t)) denoted as Constant, Step, Impulse, Gaussian, and Log-normal.
In goldfish hepatocytes, addition of ATPγ S (a non-hydrolyzable analog of ATP), activates the release of ATP following a linear function with the concentration of the analog (Pafundo et al., 2008).Assuming that ATP behaves similarly to ATPγ S, we included a factor "F " to account for a positive feedback process whereby ATPe can amplify J NL , as follows: The values of "a" and "b" are obtained by fitting a linear function to experimental data depicted in Fig. 2B.Then, J NL is modeled by the product of two independent processes, that is, a function related to the release of ATP (J R ) and the feedback mechanism (F ).
J NL applies only to viable cells, since membrane rupture gives rise to a different mechanism of ATP release as explained below.
2.1.1.2.Lytic flux Following the hypotonic challenge, a small but nevertheless significant loss of cell viability takes place.In terms of the models, cell death leads to the release of ATP to the extracellular medium by J L .This was estimated by using experimental data on the intracellular ATP content and the time course of cell viability (Fig. 2C).Cell viability, assessed by fluorescence microscopy using calcein (Pafundo et al., 2008), diminished 2.0% during the first minute of hypotonic exposure, and remained constant thereafter.The basic assumption is that the intracellular compartment contains the total mass of intracellular ATP (m ATP i ) for any given number of cells, a mass which according to the model, is released and diluted instantly into the extracellular compartment following cell death.
The ATPe concentration generated at any given time by J L can be modeled according to the following expression: Where [ATP] L e is the concentration of extracellular ATP generated by cell death (i.e., by J L ) and α is the total number of cells in the assay chamber.J L can be calculated as follows: Although this lytic flux can be considered as an important source of ATPe, it was previously demonstrated that it does not suffice to account for the experimental time course of [ATP] e (Pafundo et al., 2008).

Flux of ATP consumption by ecto-ATPase activity
We have previously identified and characterized ectonucleotidases present at the surface of goldfish hepatocytes which promote the hydrolysis of extracellular ATP, i.e., enzymes displaying ecto-ATPase activity (Schwarzbaum et al., 1998;Alleva et al., 2002;Pafundo et al., 2004).In the model, J V accounts for the flux of ATP consumption mediated by ecto-ATPase activity.The expression for J V was derived empirically after analyzing data of ectoATPase activity vs.
[ATP] e of goldfish hepatocytes.In a first approach, a hyperbolic function was fitted to experimental data reported by Pafundo et al. (2008;see 3A).
where "V m " represents the apparent maximal value of ecto-ATPase activity, and K 1/2 the concentration of extracellular ATP at which a half-maximal activity is obtained under the specific conditions of the experiment.However, since the model simulates ATPe in a concentration range with lies at least 300 times lower than the K 1/2 value, J V can be described by a linear function (Fig. 2D) as: (5)

Development of a general mathematical expression of the model 1 (without diffusion)
The model is one-dimensional and has two compartments (Fig. 1A), so that ATP can be located either in an intracellular compartment (denoted as i) or in an extracellular com-partment (denoted as e).Formally, the model can be expressed by the following ordinary differential equation: where d [ATP] e /dt represents the rate of change of [ATP] e and J NL , J L , J V are the fluxes mentioned above.
For the time course of viability (δ(t)), an exponential function was fitted to experimental data as follows (Fig. 2C): Where A 1 and A 2 are mortality and viability, both evaluated at infinitum, whereas k L is the rate coefficient of the lytic process.Thus, according to Eqs. ( 2), (3), and ( 7) J L can be fitted to a simple exponential decay function: Where A L depicts the amplitude of the lytic response calculated as . The expression for J NL includes J R , the feedback factor F and the time course of cell viability.Thus, Using Eqs. ( 1) and ( 7), we get Considering Eqs. ( 5), ( 6), (8), and (10), the following expression can be obtained: Then Eq. ( 11) can be expressed as where and Equations ( 12), ( 13), and ( 14) represent a simple mathematical model explaining the release of extracellular ATP in goldfish hepatocytes.Differential equation ( 12) governs the kinetics of [ATP] e , the state variable of the model.

Model 2 (with diffusion)
This is a three-compartment model.Similar to model 1, we have considered the ATP fluxes in one space dimension.However, unlike in model 1, here ATPe can be transported randomly by diffusion within the extracellular compartment.To build this model, we have assumed that the state variable ([ATP] e ) can follow a nonuniform distribution.This condition requires the continuous one-dimensional space to be discretized into a finite number of parts denoted as "extracellular compartments."In this respect, the simplest model would be one in which there is an intracellular compartment "i" and two extracellular compartments "e1" and "e2" (Fig. 1B).Compartment e1 corresponds to the extracellular compartment near the cell surface, whereas e2 represents the extracellular compartment surrounding e1.Then the model 2 is represented by the following equations: where [ATP] e1 and [ATP] e2 are state variables defining the extracellular concentrations of ATP in e1 and e2 compartments, respectively.In this model, [ATP] e can be calculated according to: where V 1 and V 2 are the volumes of the extracellular compartments e1 and e2, and V is the total volume resulting from the sum of V 1 and V 2 .Intracellular ATP can be released (by lytic and nonlytic means) to a small volume surrounding the cells (denoted as adjacent volume e1), with subsequent diffusion into the rest of the extracellular compartment (e2).Since EctoATPase activity is due to enzymes located on the extracellular surface of the cells, we assume that J V can only act on the adjacent volume e1, so that ).Similarly, it is supposed that the positive feedback mechanism (F ) is working on the nonlytic flux sensing the concentration of ATPe in the e1 compartment, where F = F ([ATP] e1 ) and J NL = J NL ([ATP] e1 , t).Then in model 2, J V , F , and, consequently, J NL are the same functions defined for model 1, except that [ATP] e was replaced by [ATP] e1 .Diffusion of ATPe mediates the transport of ATP between e1 and e2 according to the following discrete expression: where is the diffusive flux of ATPe, D is the diffusion coefficient of ATP in water (3 • 10 6 cm 2 sec −1 ; Hubley et al., 1996), and σ is the area of the measuring chamber (3.8 cm 2 ).1x 1 is the length of compartment e1, arbitrarily taken as the diameter of one hepatocyte (15.44 µm, Pafundo et al., 2008) and 1x 2 is the length of the most outer compartment defined as 1x 2 = V σ − 1x 1 , with V being the volume of the chamber where the experiments were performed (V = 40 µl).
By following a similar procedure as in Section 2.1.2,we arrived at a differential equation defining how [ATP] e varies in terms of the concentrations of ATPe in each compartment. and

Profiles of nonlytic release of ATP (J R )
In Eq. ( 9), we show that J NL can be obtained by multiplying the number of viable cells by J R and F (the feedback mechanism).Here, we tested different expressions (illustrated in Fig. 3) accounting for J R , so as to check which one allows the model to display the best fit to experimental data.
Types of J R functions: (a) Constant.This function is defined as Where K is a positive constant.(b) Step.Strictly, this is a Heaviside function, defined as Where K and t min are constants.The constant function can be regarded as a Step function in which t min = 0. (c) Impulse.This function is defined as Where K, t min , and t max are constants.J R has the form of a rectangular pulse that can be triggered and shut off at variable times.(d) Gaussian.
Where A G , ω G , and t G are constants.The function displays a symmetrical nonlinear increase to and decrease from a maximum value.(e) Lognormal.
Where A L , ω L , and t L are constants representing an amplitude of the flux, a time dispersion factor, and a critical time factor.This function is similar to a Gaussian expression except that the rise and decay phases might not be symmetrical.
Introduction of each of these J R functions (a to e) into general Eqs. ( 12)-( 14) as well as ( 19)-( 22) allows to test both, models 1 and 2, respectively, with five different types of J NL .It remains to be determined which is the best J NL describing the actual process of release of ATP by goldfish hepatocytes.

Model 1 (without diffusion)
In swollen goldfish hepatocytes, [ATP] e first increases to a maximum after which it decreases continuously with time toward 0. This suggests that lim t→∞ [ATP] e (t) = 0.
In order to investigate whether the model without diffusion can reproduce the experimental ATPe kinetics, we evaluated d [ATP] e /dt in Eq. ( 12) at infinite time.Since the coefficients of Eq. ( 12) are time dependent (R = R(t) and S = S(t)), and the values of R and S depend on the specific J R (t) function being used, we first evaluated each J R at infinitum.Using Eqs. ( 23) and (24) for constant and step time profiles, we find that Taking into account Eqs. ( 13) and ( 14), we get Then for time tending to infinite, the kinetics of [ATP] e will converge asymptotically to the solution of the following differential equation: Calling μ = ( V M K 1/2 100)/(aA 2 α), it is straightforward that the solution will diverge if K > μ, showing that lim t→∞ [ATP] e → ∞, a solution lacking physical meaning.In contrast, only if K < μ the solution converges at infinitum.However, even in this case, only when K = 0 the concentration of ATPe evaluated at infinitum is zero.Since in this situation J NL is null, the only source of ATPe will be J L .However, it was previously demonstrated that J L alone cannot account for the experimental results.Then model 1 loaded with Constant and Step functions cannot be used for this system.
On the other hand, Eqs. ( 25)-( 27) show that in the cases of Impulse, Gaussian, and Lognormal functions, the following behavior can be obtained at infinitum.
and using Eqs.( 13) and ( 14): Equations ( 33) and ( 34), together with Eq. ( 12) show that lim t→∞ [ATP] e = 0. Then if J R are Impulse, Gaussian, or Lognormal functions, [ATP] e will be 0 at infinitum, and consequently the model will be in a stationary state.

Model 2 (with diffusion)
As in the previous section, to analyze the [ATP] e behavior at infinite time for model 2 loaded with the Constant or Step functions the differential equations coefficients must be evaluated at infinitum, so that Eq. ( 20) can be transformed into Similarly, Eq. ( 22) turns to Calling lim t→∞ ρ(t) = β and lim t→∞ ω(t) = ε, Eq. ( 19) converges asymptotically to the following expression: Since, as mentioned in Section 3.1, experimental evidence shows that [ATP] e tends to 0 at infinite time, we investigated whether Eq. ( 37) has stationary behavior, that is, we studied Eq. ( 37) when d [ATP]e dt = 0. Thus, assuming β 6 = 0, the stationary [ATP] e in both compartments ([ATP] s e1 and [ATP] s e2 ) will be ruled by: Since [ATP] s e1 is not necessarily equal to [ATP] s e2 , the stationary [ATP] e could be nonuniformly distributed-within the extracellular compartment-at infinite time.Experimental results show that [ATP] e tends to zero at infinitum (Fig. 2A).Physically, concentrations and volumes cannot be negative.Thus, the only possibility to obtain this behavior using model 2 loaded with constant and step functions is to assume that [ATP] s e1 = [ATP] s e2 = 0, a condition that can only be achieved when K = 0. Following the reasoning shown in the previous section, the only source of ATPe in this case is the lytic flux.The fact that modeling with this single ATP flux cannot explain the experimental results makes model 2 loaded with Constant and Step functions useless.
Next, models with diffusion loaded with Impulse, Gaussian, and Lognormal function were analyzed at infinitum.Then calling lim t→∞ ρ(t) = ψ and lim t→∞ ω(t) = χ in the context of model 2 we obtained Considering these equations, expression (19) converges asymptotically to By imposing the stationary condition d[ATP]e dt = 0, and assuming that ψ 6 = 0, expression (41) gives Two conclusions can be extracted from expressions ( 38) and ( 42).First, the only possibility to obtain a stationary solution of model 2 loaded with Impulse, Gaussian, and Lognormal functions is [ATP] s e = 0, in contrast to model 2 loaded with Constant and Step functions.Secondly, since V 1 < V 2 , then θ 6 = 0. Thus, use of Impulse, Gaussian, and Lognormal J R in model 2 drives to the fact that [ATP] s e2 is zero only if [ATP] s e1 is zero, and vice versa.This case implies that [ATP] s e = 0, which is the behavior suggested by the experimental results.It is interesting to highlight that for model 2 loaded with either Impulse, Gaussian, or Lognormal functions the fact that [ATP] s e1 = [ATP] s e2 makes the [ATP] e uniform within the extracellular medium.

Preliminary information and methods
The kinetics of [ATP] e was studied by mathematical modeling.For this purpose, a computer program was developed that uses algorithms coded in visual basic to both numerically solve the differential equations involved in each type of model and fit the model to experimental data (software and source code are available).Initially, the program is loaded with one type of model including a specific J R function accounting for the nonlytic release of ATP.To perform the simulations, Eqs. ( 12)-( 14) (no diffusion) and ( 19)-( 22) (diffusion) were integrated numerically, employing the Euler method with an integration step of 1 sec.This step guaranties the order of accuracy and the stability of the Euler method, since small deviations from the true solution do not tend to grow as the solution is iterated.
The initial condition of the mathematical model was [ATP] e (t = 0) = 0.This is because in the real experiments, at time 0, goldfish hepatocytes were exposed to a hypotonic medium where ATPe was absent (Pafundo et al., 2008).
Model dependent fit to experimental data: Algebraically, the problem of fitting consists in the exploration of the free parameters space in order to minimize a given function, in this case, the sum of the squares of the residues between the experimental data of [ATP] e and that simulated by the model under study.We called the latter procedure a modeldependent fit to emphasize the fact that during the fitting the constraints imposed by the mathematical model were taken into account.In the present study, model dependent fit to experimental values of [ATP] e vs. time allowed to obtain the best values for the parameters defining J R .To make the model more versatile, the value of the ratio V m /K 1/2 (see J V in Eq. ( 5)) was allowed to vary within 1 standard error of the mean experimental value.To choose among different functions describing J R the second-order Akaike information criterion was applied (Akaike, 1992).Accordingly, after model dependent fits, the J R function that provided the lowest Akaike score (AIC) was selected.As explained in Section 2, model 1 (without diffusion) is one-dimensional and has two compartments, so that ATP can be located either in an intra or extracellular compartment (Fig. 1A).The only state variable is [ATP] e , ruled by Eqs. ( 12), ( 13), and ( 14).These equations were loaded, with each of the J R functions under study.Solving these equations results in a prediction of the time course of [ATP] e following the hypotonic insult.During the simulation, model dependent fits were applied to experimental data as follows: Five different simulations (with one type of J R function at a time) were performed where the model is fitted to experimental values of the kinetics of [ATP] e (Fig. 2A) by adjusting the values of the parameters defining J R .This procedure generates five different time profiles of J R .The model simulates [ATP] e vs. time during the hypotonic challenge as follows: dt is the concentration of ATPe at time t .Model 2 has an intracellular compartment "i" and two extracellular compartments "e1" and "e2" (Fig. 1B).The model simulates [ATP] e vs. t as follows: [ATP] e1 and [ATP] e2 are the state variables defining the concentrations of ATPe in compartments e1 and e2, while V e1 and V e2 are the volumes of the corresponding compartments, defined as V e1 = 1x e1 σ and V e2 = 1x e2 σ .

Simulation results
As predicted in Sections 3.1 and 3.2, simulations of the models assuming either a Constant release of ATP (Eq.( 23)) or a Step function (Eq.( 24)) did not converge to experimental data, regardless of the presence of one or two extracellular compartments (i.e., without or with diffusion; data not shown).
On the other hand, using an Impulse (Eq.( 25)), a Gaussian (Eq.( 26)) and a Lognormal (Eq.( 27)) functions the model reproduces qualitatively the experimental behavior.Figure 4 shows the simulated curves fitted to experimental data, with the values of best fit of the parameters shown in Table 1.
In order to choose the J R function to be used in the model, the Akaike criterion was applied (Akaike, 1992).Thus, Table 1 shows the Akaike score obtained after running the simulations with each of the J R functions in both models.Results showed that, irrespective of the model used, the Gaussian and Lognormal functions allowed a good model dependent fit, with the lognormal function providing the best option.

Analysis of model 1 (without diffusion) with lognormal flux of ATP
As mentioned above, the model that provides the best fit to experimental data is one in which there is no diffusion and where J R follows a lognormal function.
Using this lognormal model without diffusion, we analyzed the relative importance of each of the ATP fluxes that contributed to the kinetics of [ATP] e .Figure 5 shows that the lytic and nonlytic ATP fluxes displayed maximum values during the first two seconds, In all cases, the continuous black line represents the model dependent fit to experimental data, which is, in turn, represented as a gray continuous line.
Fig. 5 Time course of J NL , J L and J V predicted by the mathematical model without diffusion loaded with the lognormal time profile.J NL is plotted in a continuous black line, J V in a continuous gray line and J L in a dashed black line.
whereas ATP depletion by ecto-ATPase activity (J V ) shows a maximum after 2 minutes.Interestingly, a single, nonlinear burst of ATP is sufficient to increase [ATP] e to a maximum within the first 2 minutes.During the first 2 seconds, [ATP] e is highly governed by lytic and nonlytic ATP efflux-both potentiated by a positive feedback mechanismirrespective of ecto-ATPase activity.
The above considerations show the utility of the model in explaining the experimental results.In addition, the model can be used to predict the behavior of [ATP] e kinetics in different scenarios.Thus, to further analyze the effects of ectoATPase activity and the In all cases V m /K 1/2 ranged from 0.1 to 10 times the value of best fit.ATP max e (nM) is the maximal value of ATPe (nM), t max (sec) is the time taken to reach ATP max e and K d (sec −1 ) is the constant accounting for the decay of [ATP] e from ATP max e .
feedback loop mechanism on [ATP] e , performed new simulations using the lognormal model without diffusion with different values of "V m /K 1/2 " (Fig. 6) and with different values of "a" (Eq.( 1)) of the feedback mechanism (Fig. 7).We studied the time evolution of [ATP] e focusing on three parameters: the maximal value of ATPe (ATP max e ), the time taken to reach ATP max e (t max ) and the constant accounting for the decay of [ATP] e from ATP max e (decay constant K d ).Regarding ectoATPase activity (Fig. 6), by varying V m /K 1/2 from 0.1 to 10, we observed that K d varies directly and tmax and ATP max e inversely with V m /K 1/2 .We then studied the feedback loop by varying the best fitting value of factor "a" from 0.1 to 10 times.It can be seen that K d varies inversely, whereas both t max and ATP max e varies directly with "a" (Fig. 7).

Discussion
In goldfish hepatocytes, as in most animal cells studied so far, volume regulation in the face of osmotic gradients is an important strategy to prolong cell survival.Particularly during hypotonic exposure, cell swelling is counteracted by a process of volume downregulation termed RVD.This regulatory response has been shown to be finely regulated by extracellular ATP, which accumulates in the extracellular medium when it is released by swollen cells (Jans et al., 2002;Wang et al., 1996).Thus, to understand RVD, as well as many other biological responses (Chessell et al., 2005;Burnstock, 2006Burnstock, , 2007) ) that are mediated by endogenous ATPe, it is essential to understand the kinetics of [ATP] e .Accordingly, in this work, mathematical models were used to study how and why [ATP] e varies when goldfish hepatocytes are exposed to a hypotonic medium.The different features of the models were developed by taking into account previously published data on [ATP] e levels, ectoATPase activity, and cell viability of goldfish cells (Pafundo et al., 2008).
Experimental data had shown that when challenged by hypotonicity, goldfish hepatocytes release ATP by both lytic and nonlytic pathways.Once in the extracellular medium, ATP could further enhance the efflux of ATP, be hydrolyzed enzimatically by ectonucleotidases, and/or diffuse within the extracellular compartment.
In the present study, two models were developed which either exclude (model without diffusion) or include (model with diffusion) the diffusion of ATPe.That is, in model 1, there is no equation accounting for ATPe diffusion, since the nucleotide is assumed to diffuse instantaneously in the extracellular medium, whereas model 2 includes an explicit expression to allow ATPe diffusion.We arrived to a general mathematical form for each model (Eqs.( 12) and ( 19)).
Regarding model 2, it is interesting to analyze to what extent this model can be used to mimic the effect of an unstirred layer (USL) on the kinetics of ATPe.In a USL model adjacent to the cell membrane, there is a stagnant layer that can limit the diffusion of solutes, in our case, ATPe (Schulman and Teorell, 1938;Dainty, 1963;Dainty and House, 1966;Finkelstein, 1987).In model 2, compartment e1 would account for the USL whereas e2 stands for the bulk.It has been shown that the thickness of the USL is a dynamic parameter which depends on the diffusion coefficient of a particular solute, the viscosity of the solution, and both the osmotic and stirring fluxes (Pohl et al., 1998).Within the USL, the solute concentration is a function of the distance to the membrane (Finkelstein, 1987).
Although an assessment of the above mentioned parameters is beyond this study, we have investigated the consequences of varying the length of compartment e1 (1x 1 ), and consequently the length of e2.As we increase 1x 1 from 6 to 105 µm, the length of e2 (1x 2 ) decreases concomitantly from 99 to 0 µm.In the latter situation, model 2 becomes model 1.After running model dependent fits at different values of 1x 1 , we found no significant differences in the goodness of fit assessed by the Akaike criterion (data not shown).
These results agree well with the best model where the extracellular compartment is not fragmented, that is, model 1.Following the Occam's razor principle, since models 1 and 2 fit reasonably well to experimental data, model 1 would be better because it relies on fewer assumptions.This reasoning was quantitatively corroborated by using the Akaike criterion (see Table 1).

Cellular relevance of the functions describing ATP release (J R )
One of the main features of models 1 and 2 is the possibility to study several functions accounting for the nonlytic release of ATP.Accordingly, Eqs. ( 12) and ( 19) were loaded with functions of increasing complexity: Constant, Step, Impulse, Gaussian and Lognormal functions.
What could be the physiological correlate of these different types of ATP release?Following the hypotonic challenge, the release of ATP would be produced by activation of a putative population of ATP transport proteins that still remains unidentified in goldfish hepatocytes as well as in most cellular systems (Sabirov and Okada, 2005).
The Constant time profile would represent a hypothetical situation in which, following hypotonicity, cells sense the osmotic gradient and subsequently generate and maintain a constant ATP efflux instantaneously (from t = 0).This behavior is consistent with a synchronous and irreversible activation of channels or transporters responding to the hypotonic challenge.
The Step type function is qualitatively similar to the constant function, but a delay is imposed prior to the appearance of a constant ATP efflux.As explained by Hernández and Cristina (1998), such a time lag is to be expected between the perturbation of the extracellular osmolality, with subsequent activation of intracellular signaling pathways and the final triggering of ATP efflux.
The Impulse function is similar to the step type function, except that the release of ATP can be shut off after a period of time.This response assumes that all the putative channels or transporters are activated and shut off with equal kinetics, a simple but unlikely situation in vivo.A more realistic picture would be to use a smoother function such as a Gaussian.
This function as well as the impulse function provide symmetric time profiles, i.e., a symmetric increase and decrease of [ATP] e over time.So, to account for a possible asymmetry in the release of ATP, a lognormal expression was also studied.

Modeling the kinetics of ATPe with specific J R functions
Experimental data shows that [ATP] e , after attaining a maximum, decreases continuously until the nucleotide concentration is not significantly different from zero (Pafundo et al., 2008).
Regardless of ATPe diffusion, a model loaded with either the Constant or Step functions shows that [ATP] e evaluated at infinitum is equal to zero only when J NL = 0.However, in this case, J L would be the only source of ATPe, a situation where the experimental [ATP] e kinetics could not be quantitatively mimicked (Pafundo et al., 2008).Then the Constant and Step J R functions are not useful components of the nonlytic flux.In general, any J R whose stationary value 6 = 0 will behave as a Constant or Step function at infinite time and, therefore, will be useless when contrasting the model with experimental results.
In fact, simulations showed that irrespective of ATPe diffusion, the use of Constant and Step time profiles in the models failed to reproduce the experimental [ATP] e kinetics.
In contrast to the previous two J Rs considered, Impulse, Gaussian, and Lognormal functions yielded a null nonlytic flux of ATP at infinite time (lim t→∞ J NL = 0) always.In other words, irrespective of the parameters values of the three functions studied, these J R functions must obey lim t→∞ J R = 0 (Eq.( 33)), so that [ATP] e evaluated ad infinitum is 0, a result in agreement with the observed experimental behavior (Pafundo et al., 2008) (see Sections 3.1 and 3.2).Moreover, having discarded the Step and Constant functions, we found the use of either the Impulse, Gaussian, or Lognormal functions allowed the model to fit reasonably well to experimental data (Fig. 4).A comparison of the goodness of fit (by the Akaike criterion) showed that the use of the lognormal function in the model (compared to Gaussian and Impulse) provided the best fit to experimental values, irrespective of ATPe diffusion.
In summary, the model analysis at infinite time, together with the simulations of the model using the different J R functions, suggest that the experimentally observed kinetics of [ATP] e is compatible with a nonmonotonic nonlytic ATP release tending to [ATP] e = 0 at infinitum.This release, in turn, is compatible with a lognormal kinetic where diffusion on the extracellular media is rapid enough to be negligible during the modeling.

Model predictions
As mentioned above, the model without diffusion loaded with the lognormal J R provides the best explanation of the experimental results.We found interesting to analyze the predictive nature of this "lognormal" model under different conditions mimicking putative in vivo scenarios for goldfish hepatocytes as well as for other cell systems.Accordingly, simulations were performed considering different ectoATPase activities as well as potencies of the positive feedback mechanism.We focused on three key parameters accounting for the ATPe kinetics, i.e., the maximal value of [ATP] e (ATP max e ), the time taken to reach ATP max e (t max ) and the constant accounting for the decay of [ATP] e from ATP max e (the decay constant K d ).We first evaluated how changes in ectoATPase activity (J V of Eq. ( 5)) would affect the simulated changes of [ATP] e with time.The lognormal model predicts that a 100-fold increase of V m /K 1/2 (of J V ) would decrease both ATP max e by 31% and t max by 15-fold, whereas the decay constant (K d ) increases 10,000 times (Fig. 6).Thus, ecto-ATPase activity affects mainly the decrease and increase phases of the curve, with a relatively weaker action on the maximal [ATP] e attained.
In vivo, a change in the value of V m /K 1/2 could be obtained by changing the density of active ectoenzymes with similar kinetic properties (e.g., a cell system with partial inhibition of ecto-ATPase activity, downregulation or upregulation of functional ecto-ATPases), by cellular production of enzyme isoforms with different ATPe affinity or a combination of both.The predicted strong decrease of K d at relatively low V m /K 1/2 implies that significant ATPe concentrations would be present during longer times.In many cell systems, this could provoke desensitization of P2 receptors (Burnstock, 2007) with concomitant inactivation of the downstream processes leading to RVD or other biological responses (Chessell et al., 2005;Jans et al., 2002;Okada et al., 2001).
Another hypothetical situation can be envisaged when the potency of the positive feedback mechanism is altered, where-as explained before-J NL is a linear function of [ATP] e .By testing different values for the slope (a) of this function, the model predicts a strong increase of both ATP max e and t max , whereas the decay constant decreases slightly (Fig. 7).While the molecular mechanism enabling the feedback has not been identified in any cell system, though there are several systems where such an activating effect have been described (Anderson et al., 2004;Suadicani et al., 2006), it is interesting to speculate that a "high slope," i.e., very sensitive-feedback system would allow a small amount of ATP at the cell surface (e.g., produced by the death of a minor proportion of cells or else by diffusion of ATPe from another cell's surface) to trigger the RVD both autocrinally and paracrinally.

Concluding remarks
The models presented here constitute an attempt to understand the mathematical rules behind the non-linear, time dependent changes of [ATP] e of goldfish hepatocytes challenged by a decrease of extracellular osmolarity.
Model 1, condensed in Eq. ( 12), is theoretically simple and when used together with specific functions describing the nonlytic release of ATP, shows good agreement with experimental data.Many different physiological stimuli are able to trigger the release of ATP to the extracellular compartment.Moreover, ATPe has been shown in most cell systems studied so far to interact with specific surface receptors and be hydrolyzed by specific ectoenzymes.
Further studies could use this theoretical background as a starting point to elucidate the complex kinetics of ATP fluxes governing the kinetics of [ATP] e , the latter an important trigger of numerous biological responses.

Fig. 1
Fig. 1 Schematic summary of the mathematical models used.(A) Model without diffusion.(B) Model with diffusion.F , J NL , J L , J V , and J D stand for the positive feedback mechanism, nonlytic, lytic, Ec-toATPase activity and diffusive flux, respectively.

Fig. 2
Fig. 2 Experimental data of viability, EctoATPase activity, and ATPe of hypotonically exposed goldfish hepatocytes.All data was extracted from Pafundo et al. (2008) used with permission.(A) Time course of extracellular ATP.Levels of endogenous extracellular ATP (ATPe) of hypotonically exposed goldfish hepatocytes.In this experiment, 10 4 goldfish hepatocytes were exposed to 40 µl of hypotonic medium.Results are means ± SEM of 10 independent preparations.(B) Positive Feedback.ATP max e vs. ATPγ S concentration.Each point represents the maximal extracellular ATP found after exposing 10 4 goldfish hepatocytes to 40 µl of hypotonic medium with 0-5 µM ATPγ S. The continuous line represents the best fit by linear regression.Results are means ± SEM, n = 5-6.(C) Viability.Time course of goldfish hepatocytes viability (expressed as percentage).Results are means ± SEM, n = 60.Cells were loaded with calcein and the retention of the intracellular fluorophore was assessed in hypotonic medium each second during 45 min.Loss of fluorophore was interpreted as cell death.(D) Ecto ATPase activity.Ecto ATPase activity vs. ATP concentration (50-1000 nM).Each point represents the initial rate of ATP extinction.The line represents the best fit by linear regression.The slope of this line is the relation V m /K 1/2 used all along this work.Results are expressed as means ± SEM, n = 4.

Fig. 3
Fig.3Scheme of the time profiles of ATP release (J R ) from goldfish hepatocytes.Namely, Constant (dotted black line), Step (continuous gray line), Impulse (dashed black line), Gaussian (dotted-dashed gray line), and Lognormal (continuous black line).

Fig. 4
Fig. 4 Results of simulations using the mathematical model here developed with (B, D, F) or without (A, C, E) diffusion.(A and B) Impulse function; (C and D) Gaussian function; (E and F) Lognormal function.In all cases, the continuous black line represents the model dependent fit to experimental data, which is, in turn, represented as a gray continuous line.

Fig. 6
Fig. 6 Parameters of simulated Ecto-ATPase activity (J V ).(A) Concentration of extracellular ATP (ATPe, in nM) as a function of both time (min) and V m /K 1/2 (in arbitrary units).(B) ATP max e

Fig. 7
Fig. 7 Parameters of the feedback loop mechanism.(A) Concentration of ATPe as a function of both time (min) and a (0.1, 0.5, 1, 2, 3, 4 and 5 times the value of best fit).(B) ATP max e vs. a.(C) t max vs. a.(D) K d vs. a.In B, C, and D, a ranged from 0.1 to 10 times the value of best fit.Values of "a" have arbitrary units (A.U.).The meaning and units of ATP max e , t max and K d are given in Fig. 6.

Table 1
Parameters of the models.Results are values of best fit obtained by fitting the models to experimental data (see Fig.4).Models in the absence (model 1) and presence (model 2) of ATPe diffusion were employed, using different types of J R functions