Spatial stochastic modelling of the phosphoenolpyruvate-dependent phosphotransferase (PTS) pathway in Escherichia coli

Abstract

Motivation: Many biochemical networks involve reactions localized on the cell membrane. This can give rise to spatial gradients of the concentration of cytosolic species. Moreover, the number of membrane molecules can be small and stochastic effects can become relevant. Pathways usually consist of a complex interaction network and are characterized by a large set of parameters. The inclusion of spatial and stochastic effects is a major challenge in developing quantitative and dynamic models of pathways.

Results: We have developed a particle-based spatial stochastic method (GMP) to simulate biochemical networks in space, including fluctuations from the diffusion of particles and reactions. Gradients emerging from membrane reactions can be resolved. As case studies for the GMP method we used a simple gene expression system and the phosphoenolpyruvate:glucose phosphotransferase system pathway.

Availability: The source code for the GMP method is available at Author Webpage

Contact:  

1 INTRODUCTION

There is much evidence that stochastic effects play a key role in several biological processes, such as cell differentiation or gene regulatory networks and signalling pathways. The fluctuations in these phenomena are mostly due to biochemical processes involving a relatively small number of molecules, and the relevance of these fluctuations depends ultimately upon the underlying network architecture. Some networks components, such as switches, oscillators, positive and negative feedbacks, or transcription gene networks may display particular stochastic effects that cannot be captured using deterministic methods. The analysis of such stochastic effects can yield valuable additional insights. For example, it can help explaining part of the observed fluctuations in wet experiments, namely the internal fluctuations due to the chemical processes [intrinsic noise as in Paulsson (2004)]. The main advantage of stochastic results is the probability distribution of the steady state, trajectory or limit cycle (Arkin et al., 1998; Doubrovinski and Howard, 2005; Srivastava et al., 2002) as opposed to just the average value obtained with deterministic methods. Although effects are known to exist from theoretical analysis, only recently several leading groups have measured intracellular molecular fluctuations in vivo (Kettling et al., 1998; Arkin et al., 1998; Becskei et al., 2005; Thattai and van Oudenaarden, 2004; Doubrovinski and Howard, 2005; Pedraza and van Oudenaarden, 2005).

Many biochemical reactions do not occur in a completely homogeneous medium. For example, metabolic and signalling pathways involve reactions between cytosolic and membrane-bound molecules, causing an inhomogeneous spatial distribution of reactants. These membrane–cytosol reactions can only be sustained by the diffusion of cytosolic species, leading to an emergence of cytosolic gradients. The reaction rate and the diffusion coefficient will determine whether such gradient is significant. It is also probable that such gradients be unlikely for small prokaryotic cells but otherwise relevant in larger cells as shown in Francke et al. (2003).

Partly-membrane pathways can be modelled in the continuum and deterministic state space by partial differential equations (PDE). These models treat membrane reactions as boundary conditions on infinitely thin membranes. Alternatively a more ‘natural’ approach, adopted in particle-based models, is to account for the position in space of each particle (or a cluster of particles). Membrane-bound particles are then restricted to the membrane surface, or an adjacent sub-volume, and treated equally as cytosolic molecules with restricted mobility into the membrane subvolume. We base our work on this particular idea to tackle partly-membrane pathways.

Several models have been proposed to incorporate stochastic and spatial effects in biochemical networks, as reviewed in Takahashi et al. (2005). Currently there is no single model capable of efficiently coping with the broad range of spatial, temporal and concentration scales commonly found in biochemical networks. For instance, the highly accurate microscopic methods are too computationally demanding to simulate full pathways, while macroscopic ODE models cannot cope with spatial and stochastic phenomena. For this reason mesoscopic models may represent a plausible alternative approach and a compromise between computational efficiency and spatial and stochastic accuracy. A common approach for mesoscopic models is to tackle inhomogeneities, or spatial localization of reactions, by discretizing space into smaller homogeneous subvolumes. These subvolumes are considered to be a well-stirred system for which efficient well-stirred reaction methods can be used, either stochastic or deterministic (Weimar, 2002; Chopard et al., 1994; Stundzia and Lumsden, 1996; Baras and Malek Mansour, 1996; Ander et al., 2004; Bormann, 2001; Hattne et al., 2005).

In this article we present a method to model biochemical networks which combines the advantages of stochastic simulation and spatial discretization: the Gillespie-Multi-Particle method (GMP). Such a model can capture, within a certain range, spatial patterns and fluctuations due to diffusion and reactions. In contrast, models assuming a well-stirred medium can only capture reactional fluctuations. The GMP model is based on the discrete-in-space, operator-split and particle-based paradigms. These methods generally use dimensionless particles located in the lattice sites; diffusion takes place between lattice sites and reaction occurs locally in each site. Particularly, GMP builds on the multiparticle lattice gas automata model developed by Chopard et al. (1994). GMP employs the same diffusion process for particles, namely randomly distributing them among the nearest neighbours. For reactions, however, instead of Chopard's combinatoric method, we use the stochastic simulation algorithm (SSA) of Gillespie (1977). The separation of the reaction and diffusion processes distinguishes GMP from non-split models for the Reaction Diffusion Master Equation, where both processes are interlaced, such as the Stundzia and Lumsden (1996) approach. The core functionality of the Stundzia and Lumsden (1996) approach is a modified Gillespie method in which diffusion events for single particles are treated as a monomolecular reaction (diffusion events also follow an exponential distribution). Thus both processes are interlaced in a non-deterministic manner. In the GMP method we aim at reducing the computational cost of non-split models caused by the extra calculations of individual diffusion events. The reaction and diffusion processes are executed independently of each other (operator-split method) and the diffusion process is executed at predetermined times. Thus, the species population change in each site is due to an alternating reaction and diffusion process.

To test our method we first use a simple model of gene expression as analysed by van Zon and ten Wolde (2005). It illustrates how GMP performs in the regime where spatial fluctuations are important. Next, we apply the GMP method to a partially membrane pathway, viz., the phosphoenolpyruvate:glucose phosphotransferase system in Escherichia coli illustrated in Table 2. This case represents a real challenge for the method due to the complexity of the reactions, the spatial confinement of reactants, number of different species, and varying parameters. We validate the average results against a PDE model of the PTS pathway described in Blom and Peletier (2000) and used in Francke et al. (2003).

2 ALGORITHM

We combine the multiparticle method for diffusion (Chopard et al., 1994) and the kinetic Monte Carlo method (Gillespie, 1977) in our GMP algorithm. To tackle space we discretize the volume of the cell into a cubical lattice with L³ sites (see Fig. 1 for an illustration of the spatial discretization we use). Each lattice site holds a discrete number of dimensionless, uniformly distributed particles. Consequently, there is no need to account for the exact position of individual particles. As a result of the spatial discretization, the overall inhomogeneity of the system is limited by the coarseness of the lattice. The geometry of the membrane is represented by the lattice sites enclosing the cytosol (light boxes). These membrane sites (dark boxes) also hold cytosolic molecules, enabling membrane–cytosol reactions. To simulate the membrane as a wall, we use reflective boundary conditions. Owing to spherical symmetry, we also use pseudo-symmetric boundary conditions on the edges of the cytosol.

Fig. 1

2D slice representation of the spherical cell geometry, including the membrane sites in dark grey and a possible distribution of particles. The length of a site is λ and the neighbourhood of a cytosolic site consists of four sites. Reflective boundary conditions are used for sites on the edge of the domain.

An operator-split reaction–diffusion (RD) scheme alternates the execution of each process at predetermined intervals of time (see Algorithm 1). For each species S, its diffusion process will be executed at times $t_S=n_S{\tau }_{{\mathcal{D}}_S}$⁠, where n_S is the iteration number, and the diffusion step, τ_DS, is given by

$${\tau }_{{\mathcal{D}}_S}=\frac{1}{2d}\frac{{\lambda }^2}{D_S},$$

(1)

with d the dimensionality of the system.

Since each species may have a different diffusion coefficient, the next diffusion time is the minimum t_S that is larger than the current simulation time t_sim. Note that more than one species (even with different D_S) can occasionally have equal t_S. In this case all these species are diffused. The current simulation time t_sim is updated to t_S, and n_S is incremented by one for the diffused species.

Algorithm 1

Operator-split, Reaction–Diffusion algorithm.

The diffusion process is executed locally at every site. In this step particles are distributed randomly with equal probability among the six nearest neighbours (four are shown in Fig. 1 in black arrows). Chopard et al. showed that this diffusion rule reproduces the macroscopic diffusion equation in the limit of λ → 0, where λ is the lattice site's side length. Note that for large numbers of particles inside a lattice site it is not necessary to decide for a single particle where to diffuse, since the number of particles moving to a neighbouring site follows a normal distribution.

For the reaction process we use the SSA of Gillespie (1977), which solves the chemical master equation for well-stirred systems. This implies that the system is reaction limited at the scale of the lattice site. Note that we cannot have more spatial resolution than the lattice site size. Gillespie's method advances the system by executing the chemical reactions sequentially at reaction times τ_R. The multivariate distribution $P({\tau }_R,\mu )$1 for the time of next reaction τ_R that is of type μ is given by:

$$P({\tau }_R,\mu )=\{\begin{array}{cc}{a}_{\mu }\text{exp}(-a_0{\tau }_R\mathrm{),}\hfill & 0\le {\tau }_R<\infty ,\\ & \mu =1,\dots ,M,\hfill \\ 0\mathrm{,}\hfill & \text{otherwise}\end{array}$$

(2)

where a_μ is called propensity function, with M the total number of reaction channels and $a_0={\sum }_{\mu =1}^M{a}_{\mu }$⁠, where a_μ = k_μ,n^A,n^B/(V,N^A),k_μ is the reaction rate, $n^{A_i}$ is the number of particles in the site for species A, V is the site volume, and N_A Avogadro's number. We obtain τ_R and μ for a given species S as follows:

$${\tau }_R=1/a_0\text{ln}(1/r_1)\text{and}{\displaystyle \sum _{i=0}^{\mu -1}}{a}_i\le a_0{r}_2<{\displaystyle \sum _{i=\mu }^M}{a}_i,$$

(3)

where r₁ and r₂ are uniformly distributed random numbers from the interval [0, 1]. Then reaction μ is effectuated and the simulation time t_sim is increased with τ_R.

Because of the operator-splitting, the reaction process executes reaction only for the duration of a diffusion interval (while $\sum {\tau }_R<{\tau }_{{\mathcal{D}}_S}$⁠).

2.1 Choosing a lattice size

By adjusting the lattice size, λ, GMP's spatial resolution capabilities range from those offered by spatially resolved particle methods to well-stirred methods. However there is a limitation for the minimum lattice size λ. Despite the convergence of the diffusion process as shown by Chopard et al. (1994), the reaction–diffusion process does not converge for λ → 0. In that limit no reactions can occur as the probability of finding two particles in a site is zero. Approaching this limit also implies that the maximum possible reaction rate decreases. Thus, for a pair of reacting particles and a diffusion limited process, there is a minimal lattice size λ_o = 2σ, where σ is the diameter of a particle. Note, however, that under this condition GMP would act almost as Brownian dynamics with reaction. For complex biochemical systems, in general we do not find purely diffusion limited reactions, which allows us to take larger subvolumes that are considered to be well-stirred. A reasonable size of the well-stirred subvolume is given by λ_rmfp, the reaction mean free path2, which is the average distance traveled between two reacting collisions (cf. also the discussion on the validity of the reaction diffusion master equation in Baras and Malek Mansour (1996)).

3 CASE STUDIES

We validate GMP against two problems that highlight two different aspects of GMP's capabilities. First, a generic gene expression model to illustrate the effects of the lattice size on the variance of the fluctuations. Second, the PTS pathway, which includes spatially confined membrane-bound reactions resulting in the inhomogeneity of the system due to cytosolic gradients.

3.1 Gene expression

The generic model for gene expression, as analysed by van Zon and ten Wolde (2005), enables us to study noise as caused by a low number of molecules and diffusion limited reactions. In van Zon and ten Wolde detailed continuous-in-space simulations are provided of the system which we use to validate our GMP method and check our reasoning to choose a proper lattice size.

The system we look at is a closed volume V of 1 µm³ with a DNA binding site fixed in the centre surrounded by freely diffusing RNA (RNAp) polymerase. Once the DNA-RNAp complex is formed with association rate k_a, it can either dissociate back to separate DNA and polymerase (with rate k_d) or a protein P can be produced with a production rate k_prod with subsequent complex dissociation. The protein can further decay at the rate k_dec. Obviously the single protein production step in this model encompasses both transcription and translation which in fact consist of many biochemical reactions. Protein degradation is also simplified and treated as a first-order reaction. For the chemical reaction scheme and the parameters we refer to Table 1.

Table 1

Reaction scheme and parameters associated with the gene expression model

					Initial
	Reaction		Rate	Species	Concentration	Molecules	Dif.coef
DNA+RNAp	$\stackrel{k_{\text{a}}}{⟶}$	DNA-RNAp	3 × 109 M−1s−1	DNA		1	0
DNA+RNAp	$\stackrel{k_{\text{d}}}{\leftarrow }$	DNA-RNAp	21.5 s−1	RNAp	30 nM	18	1 µm2 s−1
DNA-RNAp	$\stackrel{k_{\text{prod}}}{⟶}$	P+DNA+RNAp	89.55 s−1
P	$\stackrel{k_{\text{dcc}}}{\leftarrow }$	∅	0.04 s–1

For choosing our lattice size we use the reasoning developed in Section 2.1. In van Zon and ten Wolde (2005) the size of the particles is 5 nm. We therefore obtain for the minimum λ a value of λ_o = 10 nm with L = 100. The λ_rmfp corresponds with L ≈ 18.

3.2 The PTS pathway and parameter setting

To test our GMP method in a more complex biological scenario, we have used the phosphoenolpyruvate:glucose phosphotransferase system (PTS system), as illustrated in the figure in Table 2. The PTS was chosen because of its moderate network complexity, involving 17 species, fast and slow reactions, and especially because of the membrane-cytosol reactions. The PTS is a group-transfer pathway situated at the beginning of the glycolysis pathway. It is present in many bacteria and it is involved in the uptake and concomitant phosphorylation of external glucose that feeds the glycolysis. A phosphoryl group derived from PEP is carried along a series of bimolecular reactions to the membrane-bound IICB, which is responsible for taking the external glucose and transferring it into the cytoplasm. Table 2 summarizes the complete sequence of reactions as in Francke et al. (2003). In that paper the kinetic parameters were obtained from Rohwer et al. (2000) and the diffusion coefficients were based on measurements by Elowitz et al. (1999).

Table 2

Reaction list and parameters associated to the PTS model illustrated below

Bulk concentrations have been converted into numbers of molecules using one-eighth of a sphere's volume of 1 µm radius, discretized in a lattice of L = 20. Surface concentrations of membrane-bound molecules require special consideration. We assume that the concentration reported in Rohwer et al. (2000) for IICB was calculated as if the membrane-bound molecules were diluted in the cytosol and its concentration, 15 µmM, was measured. According to this, the number of IICB molecules should be 3100. Since the procedure of fitting IICB bulk concentration in the original model by Rohwer does not clearly indicate the number of membrane molecules, we have adjusted the number of IICB molecules in such a way to approximate the steady state of IIA · P · IICB and IICB · P · Glc. In addition, since we have to compare surface concentrations and numbers of molecules, we scale the surface concentrations of the PDE solution in such a way that the initial IICB concentration matches the number of particles, as is used in Fig. 2(c). We found that 2100 particles for IICB, corresponding to an average of 4.3 particles per site, yield a good approximation for IIA · P · IICB and IICB · P · Glc but a less accurate approximation of IICB and IICB · P. Finally, external Glc molecules are confined to membrane sites in the reported concentration.

For the PTS system at steady state, the average value for λ_rmfp for cytosolic reactions is 0.06 ± 0.05. Since all species must share the same lattice size, L = 20 is a reasonable size.

4 RESULTS

4.1 Gene expression

With the parameters given in Table 1, simulations of the simple gene expression model yield an average level of protein P of 1032 ± 64, which is consistent with GFRD. However, the noise η in the levels of protein expression varies (see Table 3). We know that GMP converges to Gillespie for L = 1, and we observe the noise η converging to the value given by Gillespie. Using a lattice size which corresponds approximately with the reaction mean free path (L = 18), we obtain a noise level close to the GFRD reference value given by van Zon and ten Wolde (2005) using GFRD. By increasing the number of sites we obtain a better approximation, a noise level within 15% of the GFRD value, at the cost of more computations. However, L cannot be increased indefinitely since we then artificially reduce the probability of particle encounters and thus reduce the frequency of reactive collisions.

4.2 PTS Pathway

We have applied the GMP method to model the phosphoenolpyruvate:glucose phosphotransferase system (PTS) system as indicated in section 3.2. A first check of the validity of our method consists of looking at the System's Mass Evolution in Time (SMET) which shows the total number of particles present in the system over time regardless of their position. In Figure 2 we show the SMET of molecules, showing a good agreement with the reference solution of the PDE model by Blom and Peletier (2000). These results are a direct consequence of the initial decision to adjust the number of molecules of IICB to match the IIA · P · IICB steady state. Nevertheless, we note in Figure 2a, a very good agreement of the EI-related molecules. In Figure 2b, despite a good initial agreement (t < 0.005 min), the total number of molecules starts to diverge, and settling at levels ∼6% for HPr species, and 20% for IIA · P. Finally, the level of the phosphorylated IICB · P, is four times as high as the PDE solution. However, in absolute numbers, the extra ≈ 100 molecules in IICB · P, seem to be lacking in IICB.

Fig. 2

System Mass Evolution in Time (SMET) graph for the PTS pathway with parameters shown in Table 2 and a cell's radius of 1 µm. Only one trial is shown here for the GMP result. EI-related enzymes show a very good agreement between the GMP result and the PDE solution. The membrane concentration of IICB enzyme in subfigure (c) was scaled to match the number of molecules at t $=$ 1e-5.

In Figure 3 we show the radial gradients, at steady state (t_sim = 0.005 min) for some enzymes. To obtain the gradient along the radius we average the number of particles over all sites with the centre in the range [r − 0.5/L, r + 0.5/L). We then average again over 50 trials of the same systems with uniform random distributed particles. Therefore, results near the membrane will always contain far more points than those near the centre, so the average near the centre may show a larger deviation from the real mean. In Figure 3a, we observe two different behaviours of the obtained gradients: the phosphorylated enzymes HPr show only a difference in the gradient higher near the membrane (10%), whereas for IIA · P its gradient is upward shifted along the whole radius (30% near the membrane and 10% near the centre of the cell). These differences are linked with those observed in the SMET in Figure 2. Note, that owing to the spherical geometry, a slight divergence near the membrane has a larger contribution to the SMET difference than a divergence near the centre. Turning to the three basic enzymes shown in Figure 3b, we find for HPr a relative concentration difference between the centre and the membrane of the cell of 66% in both the PDE solution and the GMP averaged solution. However, in absolute value, GMP's solution for the radial gradient is consistently 16% lower than the PDE solution. Both the gradients for IIA and EI show an excellent agreement. Nevertheless, some of the molecules involved in their reactions, such as IIA · P, IICB · P and HPr · P, show a noticeable difference as seen in Figure 2.

Fig. 3

Profile gradients along the radius of a 1/8th of an sphere for the PTS system. Lines without glyphs are the solution of the PDE model in Blom and Peletier (2000). Each point is the average of all sites with centre in [r − 0.5/L, r$+$ 0.5/L], and then averaged again over 50 trials. Note the good agreement of IIA, a consequence of the tuning of IICB molecules to match IIA · P · IICB as seen in Figure 2c. EI has no significant gradient, while HPr differs by ∼15% from the value predicted by the PDE model. The first point near the centre shows a large deviation because it is only the average of one site and 50 trials.

5 DISCUSSION

We have developed the GMP reaction–diffusion method along the lines of other spatial stochastic kinetic model solvers, such as Stundzia and Lumsden (1996), Ander et al. (2004), Bormann (2001), Hattne et al., (2005). For GMP, the distinguishing feature is the separation of the reaction and diffusion processes at the scale of a lattice site. The consequent locality of the operations results in a fast algorithm. Since at the lattice site scale no spatial information of particles is used, GMP is very suitable for parallel computing. Such computational power will be required when simulation studies of whole-cell processes become necessary as experimental biology advances. GMP can cope with both reaction limited and diffusion limited problems, as well as with confined distributions of chemicals in biological systems. For complex biochemical networks a trade-off in choosing the lattice size is necessary, since different reaction rates and diffusion coefficients result in different optimal lattice sizes. Consequently, some accuracy, compared with (hard-sphere) particle models, is lost in order to gain computational speed. Nevertheless, the results can be used to guide experimentalists since such loss can be controlled and minimized.

The GMP method takes into account the principal sources of fluctuations, but the class of lattice-discrete methods, to which GMP belongs, smoothes out some noise if the lattice size is too large, because of the use of Gillespie as a reaction method (cf. Table 3). Note, however, that λ_rmfp is a safe lattice size, which is valid for small concentrations. For larger concentrations λ_rmfp decreases and in that case the slope of the gradients and the error we allow determine the optimal lattice size. For more complex biochemical systems there is no single optimal lattice size. In those cases an average λ_rmfp can be used as an orientative value, as we did for the PTS case in Section 4.2.

Table 3

Measurements of noise for the gene expression using GMP with various lattice sizes L

L	Gillespie	5	10	20	30	50	GFRD
η	0.028	0.047	0.039	0.11	0.12	0.17	0.14

Initial conditions and parameters are given in Table 1. The optimal L corresponding to the reaction mean free path is 18.

For the PTS system, the main focus was on the reactions between membrane and cytosolic species, as well as on the gradients that may arise due to the spatial localization of the membrane-bound reactions and the diffusion of cytosolic molecules. Owing to the relatively high concentrations fluctuations are not significant. As shown in Figure 3, only the species related to HPr and IIA show a significant gradient, in contrast to a flat gradient by EI which has the lowest concentration. This case illustrates that spatial inhomogeneities do not depend on the concentration levels, but on the whole pathway architecture, consisting of reaction rates, diffusion coefficients, spatial restrictions and dependencies.

Having tuned the initial number of membrane-bound molecules of IICB, we obtained a good approximation of IIA · P · IICB and IIA, but IICB, IICB · P and IIA · P differ from the PDE solution (see Figs 2 and 3). This difference is propagated down the pathway irregularly. For instance, it relatively strongly affects the gradient of HPr, though not in absolute sense (< 2e-7μM) (see Fig. 3b) and not of IIA. The reaction procedure between bulk (cytosol) and surface (membrane) molecules, which is approximated by a bulk to bulk reaction in a small volume, remains a topic for further research. Note, that in the GMP method we have to use a constant number of particles rather than a constant concentration for the membrane-bound molecules as used in PDE models. Although this is a more realistic representation of the underlying biological process and avoids the arbitrary use of concentrations of membrane-bound molecules, this introduces a practical difficulty, since this number will usually be unknown. In PDE models this concentration is estimated from the bulk concentration, which does not necessarily provide a correct estimation.

The GMP method can be easily extended to include diffusion over the membrane of membrane-bound particles and the effect of clusters of membrane-bound particles on a biochemical network can be studied. Other future aims are to improve the geometrical representation of membranes and the particle-based representation of cytosol-membrane reactions, and to extend the analysis of the choice of a proper lattice size.

This research was supported by the Netherlands Organisation for Scientific Research project NWO-CLS-635.100.007 as part of the Computational Life Sciences project.

Conflict of Interest: none declared.

¹  P(τ, μ)dτ is the probability that, given the state X_S at time t, the next reaction μ will occur in the infinitesimal time interval (t + τ, t + τ + dτ).

²The reaction mean free path is defined as λ_rmfp = (〈τ_R〉D_AB)^1/2, where 〈τ_R〉 is the mean reaction time given by 1/a_μ, and D_AB = D_A + D_B is the mutual diffusion coefficient.

Author notes

Associate Editor: Alvis Brazma