Introduction

Concentration gradients inside the cytosol are a vital piece of the cell's machinery. They underlie morphogenesis [1], and are involved in cell migration [2], cell growth [3], cell division [4], and the mechanical responses that follow fertilization, among other functions [5]. Stationary concentration gradients can spontaneously emerge from a homogeneous background in the presence of instabilities. In 1952 Alan Turing proposed a hypothetical sequential route for cell differentiation, based on a reaction-diffusion process [6]. Starting from a spatially uniform state he proved that stable inhomogeneities in “morphogen” concentration could spontaneously emerge through a diffusion-driven symmetry-breaking instability. His highly idealized model stimulated a large amount of work which eventually led to the observation of Turing patterns in open chemical reactors [7], [8]. Nevertheless, the occurrence of Turing patterns has not been unequivocally proven for a biochemical reaction in which reaction rates, concentrations and diffusion coefficients are within realistic physiological values. We do so in this paper. More specifically, we show the existence of cell-sized Turing patterns in a model of glycolysis, using realistic parameter values and equal diffusion coefficients of ATP and ADP.

The general concept behind Turing patterns involves a combination of short-range activation and large-range inhibition [9]. Several reaction-diffusion models have been proposed in the literature to explain, among others, coat patterns in mammals [9], skin patterns in fish [10] or shell patterns in mollusks [11]. However, none of these models identifies the “morphogens” involved and some of the reaction schemes lack a biological motivation. This gives almost unrestricted freedom to choose kinetic parameters and diffusion ratios. Thus, the astonishing visual similarities often observed between the predicted and the real patterns might be unrelated to the actual biological phenomena.

One of the early attempts to look for an example of a symmetry breaking mechanism in biology was due to I. Prigogine et al. [12]. They analyzed a spatially extended version of the (adiabatically reduced) 2-variable Selkov model of glycolysis, derived in [13] (see also [14], [15]). Namely, they added diffusion terms to each of the two differential equations describing the time evolution of the concentrations of ATP and ADP. It followed from their study that Turing patterns could only exist provided that ATP and ADP diffused at sufficiently unequal rates. This condition seems difficult to be met a priori, given the similarities between ATP and ADP. More recently, Hasslacher et al. [16] presented numerical simulations of Turing pattern formation in a closely related model [17] that were obtained assuming that ATP diffused 25 times faster than ADP. The authors gave a qualitative justification for their choice of diffusivities making an analogy between the role of enzymes in the glycolytic pathway and that of the immobile color indicators used for visualization in open chemical reactors [7], [8] which effectively reduce the transport rate of one of the species involved in the reactions [18]. However, the argument that applied to the latter case cannot be translated to the case of the glycolytic pathway without a deeper analysis. Firstly, enzymes affect both the dynamics of ATP and ADP. Thus, the rescaling of diffusion coefficients may occur in such a way so as to prevent the formation of patterns. Second, when the 2-variable reaction-diffusion model is obtained by the adiabatic reduction of the larger set of equations in which the enzyme dynamics is explicitly considered [19], the elimination of the fast variables (the enzymes) introduces changes of the same order of magnitude in both the diffusion and the reaction terms [20]. Thus, both changes need to be considered simultaneously in order to analyze Turing pattern formation in this setting. Finally, this theoretical discussion remains meaningless unless reasonably sized patterns can be shown to exist for realistic values of the reaction rates and of the free diffusivities.

In this paper we overcome the drawbacks of these previous works by showing that the full 5-variable Selkov model describing the PFK-controlled steps of the glycolytic pathway supports Turing patterns for realistic reaction rates and equal diffusivities of ATP and ADP. We show that Turing structures of subcellular size (10 µm) may be found by increasing the glycolytic flux and the enzyme concentration while keeping fixed the set of kinetic constants that give good agreement with the classic experimental results on homogeneous oscillations in yeast extracts [21], [22], [23], [24], [14]. The emergence of these patterns can be traced back to the differential interactions of ATP and ADP with PFK and its complexes. We then conclude that the key enzymatic step responsible for glycolytic oscillations may also provide a robust mechanism for the formation of steady state inhomogeneities in the concentration of ATP and ADP at the cellular and supracellular level.

Results

The 5-variable Selkov model reads [13]:(1)where E is the free enzyme (PFK), which can form the complexes, ES₂^γ and S₁ES₂^γ; the substrate, S₁ (ATP) is supplied by an external source at the rate, ν₁, and the product, S₂^γ (ADP) is removed by a first order reaction at rate ν₂[S₂]. From scheme (1) it is clear that the enzyme is inactive unless it has γ product molecules bound, forming the complex ES₂^γ. As done in [13], we set γ = 2 all throughout this paper. It is this activating step of the allosteric enzyme which provides the mechanism responsible for the instabilities. We assume that ATP and ADP diffuse while the enzyme and its complexes are immobile, due to the larger mass. We also assume that, initially, all concentrations are spatially uniform. Using mass action kinetics, the reaction-diffusion equations describing this model can be written, in dimensionless form, as:(2)where the dimensionless concentrations are:(3)with e₀, the total concentration of enzyme which remains constant and uniform throughout the evolution. The other quantities are(4)The dimensionless space and time coordinates are obtained by dividing the dimensional ones by an arbitrary length-scale, L, and(5)respectively. Assuming that the diffusion coefficients of ATP (D₁) and ADP (D₂) are equal (D≡D₁ = D₂), we take L≡10(DT)^1/2. Using the diffusion coefficients in vivo, D∼150 µm² s⁻¹ [25], we obtain L = 122 µm (T s⁻¹)^½. while the dimensionless diffusion coefficients are d₁ = d₂ = 0.01. The relatively small amount of enzyme used in experiments implies that ε≪1. For this reason, a quasi-steady state approximation of Eqs. (2) was analyzed in [13] considering, furthermore, spatially uniform concentrations.

We now constrain the values of the various parameters of Eqs. (2) according to previous estimations and measurements [13], [15], [26], [21], [22], [23]. Under the experimental conditions of [22], in which experiments are done using yeast extracts, the values of ν₁ and ν₂ at which the oscillations start are ν₁*∼5.8 µMs⁻¹ and ν₂*∼0.04 s⁻¹ [15]. The enzyme is very diluted in these experiments, with e₀* between 3 and 10 µM, the average concentrations of ATP and ADP are [S₁*] = 630 µM and [S₂*] = 150 µM, respectively, while the period of the oscillations, T, is between 3 and 5 min [15]. It has been shown unequivocally that the transition to oscillatory behavior in glycolysis is due to a Hopf bifurcation [24]. This means that the transition can be generically encountered by varying only one parameter. Equations (2) have a single homogeneous fixed point solution that depends on several parameters, which indeed undergoes a Hopf bifurcation. Thus, there are various ways by which this bifurcation can be reached. Among them, we choose the set of parameter values that are compatible with the observed frequency of oscillations and ATP and ADP concentrations. In particular, we find that at η = η* = 0.15, ν = ν*≈0.0041, ε = ε* = 10⁻⁶, α = 15, K₁ = 1500, K₃ = 1, a Hopf bifurcation occurs for the set of Eqs. (2) in the spatially homogeneous case. Using the definitions of these dimensionless quantities in terms of dimensional ones and the experimentally determined values, ν₁*∼5.8 µMs⁻¹ and ν₂*∼0.04 s⁻¹, we obtain that the total amount of enzyme at the Hopf bifurcation is e₀ = e₀*∼7.9 µM and that the various reaction rates satisfy(6)These numbers imply that, at the Hopf bifurcation, the stationary homogeneous solution of Eqs. (2) occurs at [S₁*]∼150 µM and [S₂*]∼145 µM, and that the dimensional period of the oscillations is T*∼2.7 min, which agrees with the experimentally determined values. We show the oscillatory behavior of [ATP] and [ADP] close to this bifurcation in Fig. 1 (a).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Behaviors predicted by the 5-variable Selkov model for different parameter values.

(a) Glycolytic oscillations in S₁ (solid curve) and S₂ (dashed curve) for η = 0.15, ν = 0.00345, ε = 10⁻⁶, α = 15, K₁ = 1500, K₃ = 1. (b) Linear growth rate of the unstable modes as a function of the square of the wavenumber, k for η = 1.215, ν = 0.03, ε = 0.0003, α = 15, K₁ = 1500, K₃ = 1, and d₁ = d₂ = 0.01. Inset: Evolution of [S₁] in the spatially homogeneous case for the same parameter values. (c) Turing space (shadowed domain) as a function of the (dimensionless) input and output rates of ATP (ν) and ADP (η), for the same parameter values as in (b). (d) Predicted value of the wave-length of the most unstable mode at each point in the Turing space of (c).

https://doi.org/10.1371/journal.pone.0001053.g001

Given these previous estimates, we explore the behavior of Eqs. (2), in the spatially extended case, varying the parameters ν, η, and ε and keeping the purely kinetic constants, α, K₁ and K₃ fixed at the previously mentioned values. We obtain the set of parameter values for which the homogeneous stationary solution of Eqs. (1) is unstable against spatially inhomogeneous perturbations (i.e., the Turing space) by analyzing the dispersion relation of the linearized evolution equations. We show in Fig. 1(b) the (dimensionless) linear growth rate, q, of the unstable modes as a function of the square of the (dimensionless) wavenumber, k, for ν = 0.03, η = 1.215, ε  = 0.0003. There is a bounded band of unstable modes for k≠0, while for k = 0 (see the inset) the fixed point is stable. We show in Fig. 1(c) the Turing space on the (ν,η) plane for ε = 0.0003. As it may be observed, patterns may exist for larger values of ν₁ as ν₂ also gets larger, i.e., as the glycolytic flux increases. A similar picture holds on the (ε,η) plane (data not shown): no pattern is possible at low enzyme concentration (if ε≤10⁻⁴). We show in Fig. 1(d) the wave-length of the most unstable mode, l_c, when η and ν are varied in the shaded region of Fig. 1(c). We then conclude that the characteristic size gets smaller as the rate of product removal, ν₂, becomes larger. We observe that, for the parameter values considered, the length-scale is always less than 23 µm, so that it can fit inside a typical cell. Using the definition of ε with the value of K₁ = 1500 and the previous rate constant estimates (Eqs. (6)) obtained at the Hopf bifurcation, which we assume remain fixed, we conclude that e₀≥800 µM for the Turing instability to occur. As we discuss later, the patterns arise due to the effective rescaling of the diffusion coefficients of ADP and ATP that the enzyme produces. This rescaling gets smeared out as e₀ decreases, leading, in turn, to the disappearance of the patterns.

We finally integrate numerically Eqs. (2) in a square domain of side 8L = 84.5 µm with periodic boundary conditions, using a finite-difference scheme on a 150×150 square grid. Fig. 2 shows the resulting pattern in S₁ after 10 min have passed since the initial situation, at which the concentrations are given by the homogeneous stationary solution ([S₁] = 2.4 mM and [S₂] = 71 µM), with 10% added noise. The typical size of the pattern agrees with the critical wavelength of the linear stability analysis, l_c ≈12 µm. The simulation shows that Eqs. (2) support stable Turing patterns with equal diffusion coefficients of ATP and ADP. The existence of these patterns can be understood in terms of the reduced set of equations obtained in the quasi-steady state approximation [19], [20], [27]. The effective diffusion coefficients of ATP and ADP that are obtained in this approximation are different due to the different interaction of ATP and ADP with the enzymes that belong to the pathway. In this scheme, ATP cannot bind to the immobile activated complex unless γ = 2 ADP molecules are already bound. Thus, the effective diffusivity of ADP –which plays the role of the activator of the step– is reduced by a larger amount than that of ATP.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Turing pattern obtained with the 5-variable Selkov model in two space dimensions.

Stationary pattern in [ATP] achieved after 10 min from an initial condition randomly distributed around the Turing-unstable fixed point. White corresponds to [ATP] = 2.47 mM and black to 1.1 mM. The simulation was done for η = 1.3, ν = 0.0175, ε = 0.0005, α = 15, K₁ = 1500, K₃ = 1, and d₁ = d₂ = 0.01.

https://doi.org/10.1371/journal.pone.0001053.g002

Discussion

Acknowledgments

This paper is dedicated to B. Hasslacher who recently passed away and whose pioneering ideas were a key factor for the development of the present work. We also acknowledge useful conversations with G. B. Mindlin, J.E. Pearson and A. Goldbeter.