Stochastic simulation of the spatio-temporal dynamics of reaction-diffusion systems: the case for the bicoid gradient

Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Stochastic simulation of the spatio-temporal dynamics of reaction-diffusion systems: the case for the bicoid gradient Paola Lecca1, Adaoha E. C. Ihekwaba1, Lorenzo Dematté1 and Corrado Priami1,2 1 The Microsoft Research - University of Trento, Centre for Computational and Systems Biology, Povo (Trento), Italy, http://www.cosbi.eu 2 DISI, University of Trento, Italy, http://www.disi.unitn.it Summary Reaction-diffusion systems are mathematical models that describe how the concentrations of substances distributed in space change under the influence of local chemical reactions, and diffusion which causes the substances to spread out in space. The classical representation of a reaction-diffusion system is given by semi-linear parabolic partial differential equations, whose solution predicts how diffusion causes the concentration field to change with time. This change is proportional to the diffusion coefficient. If the solute moves in a homogeneous system in thermal equilibrium, the diffusion coefficients are constants that do not depend on the local concentration of solvent and solute. However, in nonhomogeneous and structured media the assumption of constant intracellular diffusion coefficient is not necessarily valid, and, consequently, the diffusion coefficient is a function of the local concentration of solvent and solutes. In this paper we propose a stochastic model of reaction-diffusion systems, in which the diffusion coefficients are function of the local concentration, viscosity and frictional forces. We then describe the software tool Redi (REaction-DIffusion simulator) which we have developed in order to implement this model into a Gillespie-like stochastic simulation algorithm. Finally, we show the ability of our model implemented in the Redi tool to reproduce the observed gradient of the bicoid protein in the Drosophila Melanogaster embryo. With Redi, we were able to simulate with an accuracy of 1% the experimental spatio-temporal dynamics of the bicoid protein, as recorded in time-lapse experiments obtained by direct measurements of transgenic bicoidenhanced green fluorescent protein. 1 Introduction A preponderance of mesoscopic reaction-diffusion models of intracellular kinetics is usually performed on the premise that diffusion is so fast that all concentrations are maintained homogeneous in space. However, recent experimental data on intracellular diffusion constants, indicate that this supposition is not necessarily valid even for small prokaryotic cells [9, 11]. If the system is composed of a sufficiently large number of molecules, the concentration, i. e. the number of molecules per unit volume, can be represented as a continuous and differentiable variable of space and time. In this limit a reaction diffusion system can be modeled using differential equations. In an unstructured solvent, ideally behaving solutes (i.e. solutes for which solute-solute interaction are negligible) obey the Fick’s law of diffusion. However in biological systems even for purely diffusive transport phenomena classical Fickian diffusion is, at best, a first approximation [1, 2]. Spatial effects are present in many biological systems, so doi:10.2390/biecoll-jib-2010-150 1 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de that the spatially homogeneous assumption will not always hold. Examples of spatial effects include mRNA movement within the cytoplasm [12], Ash 1 mRNA localization in budding yeast [3], morphogen gradients across egg-polarity genes in the Drosophila oocyte [3], and the synapse-specificity of long-term facilitation in Aplysia [20]. The intracellular medium is not a homogeneous mixture of chemical species, but a highly structured environment partitioned into compartments in which the distribution of the biomolecules could be non-homogeneous. The description of diffusion processes in this environment has to start from a model in which diffusion coefficient contains its dependency on the local concentrations of the solutes and solvent. In order to tackle this problem, this paper presents a new model of diffusion coefficient for a non-homogeneous non-well-stirred reaction-diffusion system. In this model the diffusion coefficient explicitly depends on the local concentration of solute, frictional coefficient and temperature. In turn, the rates of diffusion of the biochemical species are expressed in terms of these concentration-dependent diffusion coefficients. In this study the purely diffusive transport phenomena of non-charged particles, and, in particular, the case in which diffusion is driven by a chemical potential gradient in the x direction only (the generalization to the three-dimensional case poses no problems) are considered. The derivation, introduced in this work, consists of five main steps: 1. calculation of the local virtual force F per molecules as the spatial derivative of the chemical potential; 2. calculation of the particles mean drift velocity in terms of F and the local frictional coefficient f ; 3. estimation of the flux J as the product of the mean drift velocity and the local concentration; 4. definition of diffusion coefficients as function of local activity and frictional coefficients and concentration; and 5. calculation of diffusion rates as the negative first spatial derivative of the flux J. Determination of the activity coefficients requires the estimation of the second virial coefficient calculated by using the Lennard-Jones potential to describe the inter-molecular interactions. The frictional coefficient is assumed to be linearly dependent on the local concentration of solute. The diffusion events are modeled as reaction events and the spatial domain of the reaction chamber is divided into cubic subvolumes of size l, that from now on will be called indifferently cells, meshes or boxes. The movement of a molecule A from box i to box j is represented by the k reaction Ai − → Aj , where Ai denotes the molecule A in the box i and Aj denotes the molecule A in the box j. The reaction-diffusion system is thus modeled as a pure reaction system in which the diffusion events are first order reactions whose rate coefficients k are expressed in terms of state-dependent diffusion coefficients. The space domain of the system is divided into a number Ns of subvolumes. The time evolution of the system is computed by a Gillespie-like algorithm [14] that at each simulation step selects in each subvolume the fastest reaction, compares the velocities of the Ns selected reactions and finally executes the reaction that is by far the fastest. To make the Gillespie approach applicable in each subvolume, the size of the mesh has to be chosen sufficiently small so that homogeneity and well-stirred assumptions on the distribution of the molecules inside are good approximations, and sufficiently large to have a number of reaction events significantly greater than one. The case study to which we applied our stochastic diffusion model is the Bicoid protein movement along the antero-posterior axis of the Drosophyla embryo. A controversy seems to be brewing over some recent theories and quantitative analyses addressing the fundamental question of how the Bicoid morphogen gradient is established and decoded in early Drosophila embryos. The numerical solution for the dynamics of morphogen diffusion from the localized doi:10.2390/biecoll-jib-2010-150 2 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de source of production and linear morphogen degradation shows that the morphogen gradient achieves a steady state over a broad spatial region on time periods one order of magnitude longer than the morphogen’s half life [4]. To demonstrate the actual participation of diffusion as a transport mechanism for Bicoid, several recent studies have evaluated the Bicoid diffusion rate [15, 16], all of which have yielded different values. As a step towards resolving this issue of disparities in data, a novel reaction-diffusion model of systems have been developed by us, where the medium in which molecules react and diffuse is highly structured, and the concentration of particles not spatially homogenous. In this model we assume that the kinetics of diffusive processes is driven by a multidimensional state- and time- dependent diffusion coefficient. We were able to simulate the dynamics of the gradient profile of the Bicoid protein in Drosophila Melanogaster embryo. Our procedure, implemented in our Reaction Diffusion (ReDi) software, reproduced with an accuracy of 1% the experimental spatio-temporal dynamics of the Bicoid, as recorded in a time-lapse experiment obtained by direct measurements of transgenic Bicoid-enhanced green fluorescent protein [16]. Using our methodology we were able to characterize the dynamic properties of the Bicoid gradient; predicting the kinetic rate of production and degradation and a range of variability of the average diffusion coefficient. We believe that these results are not only a case study to test a model and an algorithm, but they could give a contribution to the present studies aimed at defining the plausible biophysical mechanism, responsible for the formation of the Bicoid gradient; and towards the understanding of morphogenesis in developmental biology. The paper is organized as follows. Section 2 illustrates the mathematical model of the diffusion coefficients depending on the state variables of the system, and the models of virial coefficients, intrinsic viscosity and frictional coefficient are described. In Section 3 the method to estimate a suitable size for the subvolumes in which the entire reaction space has to be subdivided is explained. Section 4 describes the algorithm implementing the simulation of the model of reaction-diffusion systems. Section 5 shows the results obtained by applying the new algorithm to Bicoid protein gradient formation processes. We then end with some conclusions and future direction of the present study. 2 The model of diffusion: a generalization of Fick’s law In a chemical system the driving force for diffusion of each species is the gradient of chemical potential µ of this species. The chemical potential of any particular chemical species i is µi = µ0i + RT ln ai , (1) where µ0i is the standard chemical potential of the species i (i .e. the Gibbs energy of 1 mol of species i at a pressure of 1 bar), R = 8.314 J · K−1 · mol−1 is the ideal gas constant, and T the absolute temperature. The quantity ai is called chemical activity of component i, and it is given by γ i ci ai = 0 (2) c where γi is the activity coefficient, and c0 a reference concentration, which, for example, could be set equal to the initial concentration. The activity coefficients express the deviation of a solution from ideal thermodynamic behavior and, in general, they may depend on the concentration of all the solutes in the system. For an ideal solution, the limit of γi , which is recovered doi:10.2390/biecoll-jib-2010-150 3 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de experimentally at high dilutions is γi = 1. If the concentration of species i varies from point to point in space, then so does the chemical potential. For simplicity, here the case in which there is only a chemical potential gradient in the x direction is taken into account. Chemical potential is the free energy per mole of substance, free energy is the negative of the work, W , which a system can perform, and work is connected to force F acting on the molecules by dW = F dx. Therefore an inhomogeneous chemical potential is related to a virtual force per molecule of kB T c0 X ∂ai ∂cj 1 dµi =− Fi = − NA dx γi ci j ∂cj ∂x (3) Fdrag,i = fi vi , (4) where NA = 6.022 × 1023 mol−1 is Avogadro’s number, kB = 1.381 × 10−23 J · K−1 is the Boltzmann constant, and the sum is taken over all species in the system other than the solvent. This force is balanced by the drag force experienced by the solute (Fdrag,i ) as it moves through the solvent. Drag forces are proportional to speed. If the speed of the solute is not too high in such a way that the solvent does not exhibit turbulence, the drag force can be written as follows where fi ∝ ci is the frictional coefficient, and vi is the mean drift speed. Moreover, if the solvent is not turbulent, the flux, defined as the number of moles of solute which pass through a small surface per unit time per unit area, can be approximated as in the following Ji = ci vi , (5) i.e. the number of molecules per unit volume multiplied by the linear distance traveled per unit time. Since the virtual force on the solute is balanced by the drag force (i. e. Fdrag,i = −Fi ), the following expression for the mean drift velocity is obtained vi = Fi , fi so that Eq. (5) becomes Ji = − X kB T X ∂ai ∂cj ∂cj ≡− , Dij γi fi j ∂cj ∂x ∂x j (6) where Dij = kB T c0 ∂ai , γi fi ∂cj (7) are the diffusion coefficients. The Eq. (7) states that, in general, the flux of one species depends on the gradients of all the others, and not only on its own gradient. However, here it is supposed that the chemical activity ai depends only weakly on the concentrations of the other solutes, i.e. it is assumed that Dij ≈ 0 for i 6= j and that Fick’s laws still holds. Let Di denote Dii . It is still doi:10.2390/biecoll-jib-2010-150 4 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de generally the case that Di depends on ci in sufficiently concentrated solutions since γi (and thus ai ) has a non trivial dependence on ci [34]. It is only in one very special case, namely that of an ideal solution with γi = 1, where the diffusion coefficient, Di = kB T /fi , is constant. In order to find an analytic expression for the diffusion coefficients, Di , in terms of the concentration, ci , let us consider that the rate of change of concentration of the substance i due to diffusion is given by ∂Ji , (8) ∂x Substituting Eq. (7) into Eq. (6), and then substituting the obtained expression for Ji into Eq. (8), gives Di = − ∂ Di = − ∂x so that Di =    ∂ci − Di (ci ) , ∂x  ∂Di (ci ) ∂ci ∂ 2 ci ∂Di (ci ) ∂cj ∂ci ∂ 2 ci + Di (ci ) 2 = + Di (ci ) 2 . ∂x ∂x ∂x ∂cj ∂x ∂x ∂x Let ci,k denote the concentration of a substance i at coordinate xk , and l = xk − xk−1 the distance between adjacent mesh points. The derivative of ci with respect to x calculated at xk− 1 is 2 ci,k − ci,k−1 ∂ci ≈ . (9) ∂x xk− 1 l 2 By using Eq. (9) into Eq. (6) the diffusive flux of species i midway between the mesh points, Ji,k− 1 is obtained: 2 ci,k − ci,k−1 Ji,k− 1 = −Di,k− 1 , (10) 2 2 l where Di,k− 1 is the diffusion coefficient midway between the mesh points. 2 The rate of diffusion of substance i at mesh point k is Dik = − Ji,k+ 1 − Ji,k− 1 2 2 l , and thence Dik = Di,k− 1 2 l2 (ci,k−1 − ci,k ) − Di,k+ 1 2 l2 (ci,k+1 − ci,k ) (11) To determine completely the right-hand side of Eq. (11) is now necessary to find an expression for the activity coefficient, γi , and the frictional coefficient, fi , contained in Eq. (7) for the diffusion coefficient. In fact, by substituting Eq. (2) into Eq. (7) an expression of the diffusion coefficient in terms of activity coefficients γi is obtained Dii = doi:10.2390/biecoll-jib-2010-150 ci ∂γi  kB T  1+ fi γi ∂ci (12) 5 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Let us focus now on calculation of the activity coefficients: a way to estimate the frictional coefficients will be presented in Section 2.1. By using the subscript ’1’ to denote the solvent and ’2’ to denote the solute, it can be written that   γ 2 c2 , c0 (13) 1 1 ∂γ2  ∂µ2 . = RT + ∂c2 c2 γ2 ∂c2 (14) µ2 = µ02 + RT ln where γ2 is the activity coefficient of the solute and c2 is the concentration of the solute. Differentiating with respect to c2 gives The chemical potential of the solvent is related to the osmotic pressure (Π) by µ1 = µ01 − ΠV1 , (15) where V1 is the partial molar volume of the solvent and µ01 its standard chemical potential. Assuming V1 to be constant and differentiating µ1 with respect to c2 yield ∂Π ∂µ1 = −V1 . (16) ∂c2 ∂c2 Now, from the Gibbs-Duhem relation [31], the derivative of the chemical potential of the solute with respect to the solute concentration is M(1 − c2 v) ∂µ1 M(1 − c2 v) ∂Π ∂µ2 =− = , (17) ∂c2 V 1 c2 ∂c2 c2 ∂c2 where M is molecular mass of the solute and v is the partial molar volume of the solute divided by its molecular mass. The concentration dependence of osmotic pressure is usually written as   Π RT 2 = 1 + BMc2 + O(c2 ) . c2 M (18) where B is the second virial coefficient (see Section 2.2), and thence the derivative of Π with respect to the solute concentration is ∂Π RT = + 2RT Bc2 + O(c22 ). ∂c2 M Introducing Eq. (19) into Eq. (17) gives  1 ∂µ2 = RT (1 − c2 v) + 2BM . ∂c2 c2 From Eq. (14) and Eq. (20) it can be obtained that (19) (20) i 1 ∂γ2 1h (1 − c2 v)(1 + 2BMc2 ) − 1 , = γ2 ∂c2 c2 doi:10.2390/biecoll-jib-2010-150 6 Journal of Integrative Bioinformatics, 7(1):150, 2010 so that Z γ2′ 1 dγ2 = γ2 Z c′2 c0 http://journal.imbio.de i 1h (1 − c2 v)(1 + 2BMc2 ) − 1 dc2 . c2 On the grounds that c2 v ≪ 1 [37], solving the integral yields γ2′ = exp[2BM(c′2 − c0 )] (21) The molecular mass Mi,k of the species i in the mesh k can be expressed as the ratio between the mass, mi,k , of the species i in that mesh and the Avogadro number Mi,k = mi,k /NA . If pi is the mass of a molecule of species i and ci,k l is the number of molecules of species i in the mesh k, then the molecular mass of the solute of species i in the mesh k is given by pi l ci,k . (22) NA Substituting the expression in Eq. (21) gives, for the activity coefficient of the solute of species i in the mesh k (γi,k ), the following equation Mi,k =  pi l 2  c . γi,k = exp 2B NA i,k (23) 2.1 Intrinsic viscosity and frictional coefficient The diffusion coefficient depends on the ease with which the solute molecules can move. It is a measure of how readily a solute molecule can push aside its neighboring molecules of solvent. An important aspect of the theory of diffusion is how the magnitude of the frictional coefficient, fi , of a solute of species i and, hence, the diffusion coefficient, Di , depend on the properties of the solute and solvent molecules. Examination of well-established experimental data shows that diffusion coefficients tend to decrease as the molecular size of the solute increases. The reason is that a larger solute molecule has to push aside more solvent molecules during its progress and will therefore move slowly than a smaller molecule. A precise theory of the frictional coefficients for the diffusion phenomena in a biological context cannot be simply derived from the elementary assumptions and model of the kinetic theory of gases and liquids. Stokes’s theory considers a simple situation in which the solute molecules are so much larger than the solvent molecules that the latter can be regarded as a continuum (i. e. not having molecular character). For such a system Stokes deduced that the frictional coefficient of the (H) (H) solute molecules is fi = 6πri η, where ri is the hydrodynamical radius of the molecule and η is the viscosity of the solvent. For proteins diffusing in the cytosol, estimating the frictional coefficient through Stokes’s law is hard, for several reasons. First of all, the assumption of very large spherical molecules in a continuous solvent is not a realistic approximation for proteins moving through the cytosol: proteins may be not spherical and the solvent is not a continuum. Furthermore, in protein-protein interaction, in the cytosol, water molecules should be included explicitly, thus complicating the estimation of the hydrodynamical radius. Finally, the viscosity of the solvent, η within the cellular environment cannot be approximated either as the viscosity of liquid or the viscosity of gas. In both cases, the theory predicts a strong dependence on the temperature of the system, that has not been found in the cell, where the most significant factor in determining the frictional coefficient is the concentration of solute molecules. To model doi:10.2390/biecoll-jib-2010-150 7 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de the effects of non-ideality on the friction coefficient it is assumed that it linearly depends on the concentration of the solute, as in sedimentation processes [35]. Equation (24) gives the frictional coefficient fi,k of species i at mesh k. In this equation kf is an empirical constant, whose value can be derived from the knowledge of the ratio R = kf /[η]: fi,k = kf ci,k . (24) Accordingly to the Mark-Houwink equation [23], [η] = kM α is the intrinsic viscosity coefficient, α is related to the shape of the molecules of the solvent, and M is the molecular mass of the solute. If the molecules are spherical, the intrinsic viscosity is independent of the size of the molecules, so that α = 0. All globular proteins, regardless of their size, have essentially the same [η]. If a protein is elongated, its molecules are more effective in increasing the viscosity and [η] is larger. Values of 1.3 or higher are frequently obtained for molecules that exist in solution as extended chains. Long-chain molecules that are coiled in solution give intermediate values of α, frequently in the range from 0.6 to 0.75 [23]. For globular macromolecules, R has a value in the range of 1.4 - 1.7, with lower values for more asymmetric particles [17]. 2.2 Calculated second virial coefficient The statistical mechanics definition of the second virial coefficient is given by the following expression Z ∞ h u(r) i 2 dr (25) B = −2πNA r exp − kB T 0 where u(r), which is given in Eq. (26), is the interaction free energy between two molecules, r is the intermolecular center-center distance, kB is the Boltzman constant, and T the temperature. In this work, it is assumed that u(r) is the Lennard-Jones pair (12,6)-potential (Eq. 26), that captures the attractive nature of the Van der Waals interactions and the very short-range Born repulsion due to the overlap of the electron clouds: u(r) = 4 By expanding the term exp  1 kB T r 6 4  h 1 12 r  1 6 i . − r (26) into an infinite series, the Eq. (25) becomes Z ∞ h X 1i 1 ∗ j ∞ 2−6j (T ) r exp − T ∗ 2 dr B = −2πNA j! r 0 j=0 , where T ∗ ≡ 4/(kB T ) and thus B=− ∞  1 1  1 1 πNA X 1 j 4 (kB T )− 4 + 2 j Γ − + j 6 j=0 !j 4 2 (27) An estimate of B is given by truncating the infinite series of Γ functions to j = 4 since results not shown here prove that taking into account the additional terms, obtained for j > 4, does not significantly influence the simulation results. doi:10.2390/biecoll-jib-2010-150 8 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de 3 The optimal size of the system’s subvolumes The reaction chamber volume, V , has been divided into subvolumes of volume ∆ and side length l, on the basis of the kinetic and dynamical properties of the diffusion particles. The subvolumes have been chosen sufficiently small so that the probability distributions of the reactants can be treated as uniform inside each subvolume. This means that the rate by which two molecules in a subvolume reacts does not depend on their initial locations. Let consider diffusion as a time-dependent process, in which some distribution of concentration is established at some moment, and then allowed to disperse without replenishment. Fick’s law and its analogues for the transport of other physical properties relate to the flux under the influence of a constant gradient. They therefore describe time-independent processes and they refer, for example, to the flow of particles along a constant concentration gradient which is sustained by injecting particles in one region, and drawing them off in another. From the Fick’s second law, the mean distance through which particle of solute has spread after time t is r Dt , (28) lf = 2 π where D is the diffusion coefficient of the particle. Let te be the mean free time with respect to non-reactive (elastic) collisions and tr the mean free time with respect to reactive collisions. The net distance covered by the particle during its lifetime is s r r πlf2 tr D tr tr L=2 =2 = lf . (29) π 4te π te In order to have a homogeneous mixing inside boxes, the length l of the box side has to fulfill the following inequality, l ≪ L. (30) It is worthy of note the fact that if this inequality is fulfilled, the particles in each box obey the Einstein formula for the probability of fluctuations around the steady state. Note also that the rate at which two molecules in a subvolume react does not depend on their initial location if the inequality (30) is fulfilled. In terms of the diffusion coefficient, D, Eqs. (29) and (30) can be re-written as p l ≪ 2 Dtr . (31) Now, in order to estimate the upper bound of l the diffusion coefficient D and the reaction time tr have to be determined. The diffusion coefficient differs from species to species, and, in general, depends on the local concentration of solute. Since the local concentration of solute changes in time as consequence of the occurrence of the chemical reaction events and the diffusion events themselves, this would entail a dynamical change of l through the Eq. (31). This could make the algorithm of simulation more complex, so that, it is profitable to fix the value of l at the initialization time at p (32) l ≈ hDitr doi:10.2390/biecoll-jib-2010-150 9 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de where (dif f ) hDi = Rdif f MX 1 Di0 + Rreact i=1 (33) and Di0 is the diffusion coefficient of the i-th species at time t = 0, and M dif f is the number of species that diffuse. In the next section the model of the diffusion coefficient as function of local concentration and the waiting time of reaction, tr , are explained. 3.1 The waiting time of reaction Let Ri be the i-th reaction channel expressed as r i → ... Ri : li1 Sp(i,1) + li2 Sp(i,2) + · · · + liLi Sp(i,Li ) − where lij is the stoichiometric coefficient of reactant Sp(i,j) , p(i, j) is the index that selects the species that participates in Ri , Li is the number of reactants in Ri , and ri is the rate constant. If the fundamental hypothesis of stochastic chemical kinetics [14] holds within a box, both diffusion and reaction events waiting times are distributed according to a negative exponential distribution, so that a typical time step has size 1 tr ≈ R X R aν ν=1 −1 1 = R  RX dif f (dif f ) ai + RX react i=1 i=1 (react) ai  (34) where R is the number of events. It is given by R = Rdif f + Rreact , where Rdif f is the number of possible diffusion events and Rreact is the number of reaction events [5]. The diffusion and reaction propensities are given by the following expressions, respectively (dif f ) ai = (dif f ) ri QMi(dif f ) j=1 (#Sp(i,j) )lij QLi(dif f ) j=1 (react) ai = (react) ri QMi(react) j=1 (#Sp(i,j) )lij QLi(react) j=1 (dif f ) (react) , (35) lij ! , (36) lij ! where Mi and Mi are the number of chemical species that diffuse and the number of (dif f ) (react) those the undergo to reactions, respectively. In general M 6= Mi + Mi , since some (dif f ) species both diffuse and react. In Eq. (35), ri is the kinetic rate associated to the jumps (react) is the stochastic rate constants between neighboring subvolumes, whereas in Eq. (36), ri of the i-th reaction. From Eq. (11), the rate coefficient of the first order reaction representing a diffusion event is recognized to be as follows. (dif f ) ri doi:10.2390/biecoll-jib-2010-150 = Dii . l2 (37) 10 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de 4 The algorithm and data structure In this section we illustrate how we incorporated the state-dependent diffusion model presented in the previous section into a Gillespie-like approach of stochastic simulation of biochemical kinetics. For the reader’s convenience, we report a brief description of the Gillespie Direct and First Reaction methods [14]. Let assume that in the system there are R reactions and M chemical species. Let assume also that the solution is well mixed and in thermal equilibrium. At ~ any instant of time the system is described by the state vector X(t) = {X1 (t), . . . , XM (t)}. Gillespies algorithm addresses two questions: 1. Which reaction occurs next? 2. When does it occur? Both of these questions must be answered probabilistically by specifying the probability density P (µ, τ ) that the next reaction is µ and that it occurs at time τ . It can be shown [14] that  P (µ, τ ) = aµ exp − τ R X j=1  aj dτ. (38) This equation leads directly to the answers of the two aforementioned questions. First, the probability distribution for reactions is obtained by integrating P (µ, τ )) over all τ from 0 to ∞, i. e. aµ , (39) Pr(Reaction = µ) = PR j=1 aj aj where aj the propensity of reaction j, as in Eqs. (35) and (36). Second, the probability distribution for waiting reaction time is obtained by summing P (µ, τ ) over all τ , i. e. P (τ )dτ = R X j=1  aj aj exp  −τ R X j=1  aj dτ. (40) The next reaction and its waiting time are realizations of the reaction probability distribution, and time probability distribution, respectively. The stochastic simulation according to the Gillespie Direct method consists in the following steps. ~ 1. Set initial numbers of molecules in X(t), set t ← 0, and the absolute simulation time T . 2. Calculate the propensity functions, aµ , for all j = 1, . . . , R. 3. Choose j according to the distribution in Eq. (39). 4. Choose τ from an exponential distribution with parameter PR j=1 aj (as in Eq. (40)). 5. Change the number of molecules to reflect execution of reaction µ. Set t ← t + τ . doi:10.2390/biecoll-jib-2010-150 11 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de 6. Go to Step 2 and repeat the procedure until t ≤ T . The algorithm is direct in the sense that it generates µ and τ directly. Gillespie also developed the First Reaction Method which generates a putative time τj for each reaction to occur - a time the reaction would occur if no other reaction occurred first - then lets µ be the reaction whose putative time is first, and lets τ be the putative time τj . Formally, the algorithm for the First Reaction Method is as follows: ~ 1. Set initial numbers of molecules in X(t), set t ← 0, and the absolute simulation time T . 2. Calculate the propensity function, aµ , for all j, j = 1, . . . , R.. 3. For each µ, generate a putative time, τj , according to an exponential distribution with parameter aj . 4. Let µ be the reaction whose putative time, τj , is least. 5. Let τ be τj . 6. Change the number of molecules to reflect execution of reaction µ. Set t ← t + τ . 7. Go to Step 2 and repeat the procedure until t ≤ T . The Direct and First Reaction Methods are provably equivalent [14], that is, the probability distributions used to choose µ and τ are the same in the two methods. With regard to the complexity of the First Reaction Method, this algorithm uses R random numbers per iteration (where R is the number of reactions), takes time proportional to r to update the a’s, and takes time proportional to R to identify the smallest τj . The design of the algorithm implementing our model of reaction-diffusion system inspires to the structure of the Next sub-volume method proposed by Elf et al. [10]. This method selects the next reaction and the time at which it will occur by using the Gillespie First Reaction method [14]. Each cell and the corresponding reaction time and reaction type is stored in a global priority queue that is sorted with increasing writing reaction time. From this queue at each time step, the fastest reaction (i.e. the reaction with the smallest waiting time) is picked and executed. Once the reaction has been executed, the state of the cell, as well as the state of the neighboring cells that have been affected by the occurrence of this reaction, are updated. This approach is efficient as it does not update the state of all the cells, but only those of cells in which the occurrence of a reaction has produced changes in the inner amount of molecules. However, the method is centralized and sequential and does not scale to very large systems. Moreover, it cannot be easily adapted to turn parallel or distributed computing procedures to profit. Since the number of reactions involved in the system could be of the order of millions, the property of scalability is required to make large simulations feasible. The algorithm proposed by the authors overcomes the scalability’s limitations of the Next sub-volume method because it doe not make use a global priority queue. For each cell a set of dependency relations with neighbor cells is drawn; in a cell an event (reaction or diffusion) can be executed only if it is quicker than the diffusion events of neighbor cells, since the diffusion events in and out of the cell could change the reactant concentrations, doi:10.2390/biecoll-jib-2010-150 12 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de and, consequently the reaction propensities and the waiting times of the events in the neighbor cells. The algorithm has still the same average computational complexity of Elf’s methods. Nevertheless, by removing the global priority queue and introducing a dependency relations graph, the algorithm gains the scalability property. The new algorithm consists of the following steps. ~ 1. Set initial numbers of molecules in X(t), set t ← 0, and the absolute simulation time T . Divide the reaction chamber volume V into boxes of size l as in Eq. (32). 2. In each cell, calculate the time and the type of the next event with the FRM and store them in a private priority queue, ordered with increasing waiting time. 3. Each cell “communicates” with its neighbors in a hierarchical way, on the basis of the dependency relations, to decide which one holds the event with the smallest waiting time, say τs , that will be executed next. Execute the event and update the state of the cell and those of the neighbor cells, in the case in which the event is a diffusion. 4. Update the time variable: t ← t + τs . 5. Go to Step 2 and repeat the steps until t ≤ T . 4.1 The prototype Redi We developed the software tool Redi, in order to implement the reaction-diffusion system dynamics in the stochastic framework of the algorithm described above. This tool is part of the software platform CoSBiLab that implements a new conceptual modeling, analysis and simulation approach - primarily inspired by algorithmic systems biology [33] - to biological processes. CoSBiLab tools have been developed by our group and an overview is available at http://www.cosbi.eu/index.php/research/prototypes/overview Redi works in Windows and is freely available for non commercial purposes at the following url http://www.cosbi.eu/index.php/research/prototypes/redi It is a command line tool, that is invoked by typing on the command shell under Windows the command redi.exe [params] In Figure 1 we show an exaple of specification of the simulation parameters. The command shown in Figure 1 launches a simulation of 200000 steps, recording the state of the system every 1000 steps. The width of the simulation space is 32, the height is 32, and the thickness is 4. The length of the side of a mesh is 700 (the units depend on scale of the system under investigation). The command invokes three writers: GlobalDataWriter, which records overall system concentrations in textual form, AVFDataWriter, which records doi:10.2390/biecoll-jib-2010-150 13 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de redi.exe /steps:200000 /deltaSteps:1000 /fileName:Input_file_name.txt /outputPrefix:Redi_output /x:32 /y:32 /z:4 /cellLength:700 /writer:GlobalDataWriter /writer:AVFDataWriter /writer:Cosbi.DiffuseSim.ViewVolumeDataWriter,Simple3DView ... /writerSpecies:Species1 /writerSpecies:Species2 Figure 1: Example of command line syntax invoking Redi simulator (see the text for more explanations). concentrations of each species in each voxel in Virvo data format. The third writer is a custom data writer, provided in a plug-in included with the download package. It simple capture the output and writes it to the screen, using a 3D dithering view to show species concentrations in space. Notice that for this writer as for any other writer provided in a plug-in, it is necessary to provide a fully qualified name in the form Namespace.ClassName,AssemblyName. The assembly name is usually the name of the .dll file in which the plug-in is stored. The reference manual of Redi available on Redi web page, reports the full list of writers provided by CoSBi. It is possible to specify which chemical/biological species will be handled by the writers with the /writerSpecies switch (in the example of Figure 1 Species1 and Species2 are specified). The /fileName switch, used to pass the input model file to the simulator, is always required. Models for the reaction-diffusion system are given in a simple input language in a text input file. An example of input file is shown in Figure 2. Named biochemical entities (molecules, atoms, ions, complexes, etc) that are included in the system are declared in the first section of the input file, preceded by the keyword var (lines 3-7 in Table 2). For each entity it is possible to specify (i) a constant diffusion rate (preceded by the keyword rate), only if their diffusion is independent of the local conditions of the medium; in this case the dynamics of the diffusion reduces to Fick’s laws. (ii) The molecular mass, expressed in kDa, after the keyword weight; in this case the diffusion of species is state-dependent and its dynamics is predicted by our model of Fick’s generalized law. The algorithm takes as input the molecular mass of the species to compute its chemical activity (see Eq. (22)), and uses it to compute Dii with the Eq. (12); finally, the rate parameter of the first order reaction representing the diffusion event of the molecule moving from one cell to another is computed by Eq. (37). Nothing is specified after the declaration of species D and E, that means that D and E are not involved in the diffusion transport, so the user does not need to indicate value the diffusion coefficient or the molecular mass. Reactions can be specified in a very intuitive format (lines 12 and 13 in Figure 2): at line 12, a bimolecular reactions is specified. The rate constant is given in scientific notation between square brackets. The bimolecular reaction at line 13 is a shorthand to specify reversible reactions. In this case, it is necessary to specify between the square brackets the rate constant for the forward and then the rate constant for the backward reactions. The reaction list has to be followed by the specification of the concentration (or number of molecules) and spatial location of each species on the grid (lines 21-37). After the species name, doi:10.2390/biecoll-jib-2010-150 14 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de the first three numbers within the square brackets indicate the spatial coordinates (x, y, z), and the fourth number is the abundance of the species in (x, y, z). Although this input format may be verbose, it allows for great control and precision. We are currently studying better ways of inputting both data (location and amount of particles) and geometry (that is, for now, fixed and with a regular polyhedral shape). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 // Species in the systems var var var var var A : rate 10e-3; B : rate 10e-6; C : weight 24.0; D; E; // Reactions B + C -> A + C [1.9E-03]; D + B -> E [5.9E-03, 20]; run // Initial conditions: // species amount in the cells of the spatial grid A [253, 33, 0, A [256, 33, 0, A [262, 33, 0, ... B [257, 33, 0, B [260, 33, 0, ... C [245, 35, 0, C [246, 35, 0, C [247, 35, 0, ... D [245, 35, 0, D [346, 35, 0, ... E [253, 33, 0, E [262, 33, 0, ... 17]; 0]; 7]; 2]; 25]; 47]; 3]; 36]; 13]; 13]; 0]; 9]; Figure 2: Example of an input file for the Redi simulator. Starting from the initial conditions, the stochastic algorithm of Redi selects the event (reaction or diffusion) that has to occur accordingly to its waiting time, executes it and updates the current doi:10.2390/biecoll-jib-2010-150 15 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de state of the system defined by the abundance of chemical/biological species in each sub-volume of the space. The Redi package includes a number of modules providing different kinds of output, e. g. text files, writing data files to be used with professional visualization packages, visualizing 3D concentration gradients, etc.) as well as a simple interface to write custom output modules (see some example of output visualization in Figs. 3 - 5). Figure 3: Redi’s screenshots of output’s visualization. Figure 4: Visualization of Redi outputs: x − y scatter-plot of the time behavior of average concentrations of species in the system (left lower corner) and 3D visualization of the concentration gradient at a certain step of simulation (on the right). The menu on the right upper corner lists the parameters of the simulation that are tunable by the user, such as the number of sub-volumes, their size, the number of steps of the simulation etc.. 5 Case study: modeling the formation of Bicoid gradient One of the most fundamental problems in developmental biology is the question of morphogenesis; that is, by what mechanisms do cells self-organize into highly complex spatial distributions doi:10.2390/biecoll-jib-2010-150 16 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Figure 5: Visualization of Redi outputs as heatmaps. The Red color indicates regions with a low concentration of molecules, while yellow and white colors indicate regions occupied by large amounts of molecule concetration. giving rise to tissues and organs during embryonic development? While tremendous advances have been made in the last fifty years, there are still unresolved questions on a morphogen’s mechanism of transport and relationships. The first and perhaps still most elegant example of a morphogen is bicoid. The gradient formation of which acts as a suitable case study of a reaction-diffusion system; where the medium, that is the syncytial embryo, is inhomogenous and highly dynamic, and the most pronounced changes associated with the number and the spatial distribution of the nuclei [7]. However, there seems to be brewing a controversy over some recent theories and quantitative analyses addressing the fundamental question of how the bicoid morphogen gradient is established and decoded in early Drosophila embryos. The numerical solution for the dynamics of morphogen diffusion from a localized source of production and linear morphogen degradation shows that the morphogen gradient achieves a steady state over a broad spatial region on time periods one order of magnitude larger than the morphogen’s half life [4]. To demonstrate the actual participation of diffusion as a transport mechanism for the bicoid protein, several recent studies have evaluated the bicoid’s diffusion rate (see the works of T. Gregor et al. in [15, 16] for a comprehensive review of the literature on this debate), all of which have yielded different values. As a step towards resolving this issue of disparities in data, a novel reaction-diffusion model have been developed by us, where the medium in which molecules react and diffuse is highly structured, and the concentration of particles not spatially homogenous. In this model we assume that the kinetics of diffusive processes are driven by multidimensional state- and time- dependent diffusion coefficients. We were able to simulate the dynamics of the gradient profile of the bicoid protein in the Drosophila Melanogaster embryo. Our procedure, implemented in our Redi software, reproduced with an accuracy of 1% the experimental spatio-temporal dynamics of bicoid, as recorded in a time-lapse experiment obtained by direct measurements of transgenic bicoid-enhanced green fluorescent protein [16]. doi:10.2390/biecoll-jib-2010-150 17 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Using our reaction-diffusion model we were able to characterize the dynamic properties of the bicoid gradient; predicting the kinetic rate of production and degradation and the range of variability of the average diffusion coefficient. We believe that these results could give a contribution to the present studies aimed at defining the plausible biophysical mechanism responsible for the formation of the bicoid protein gradient, and towards the understanding of morphogenesis in developmental biology. 5.1 The Redi model of the Bicoid gradient In Drosophila Melanogaster, bicoid is a maternally transcribed gene that organizes anterior development. Its mRNA is localized at the anterior pole of the oocyte and is translated soon after fertilization. As a consequence of the anterior localization of mRNA a gradient of the bicoid protein is formed along the antero-posterior axis, simultaneously with nuclear cleavage cycles. Because the bicoid protein is a transcription factor that functions when the early Drosophila melanogaster embryo is a syncytium, the bicoid protein confers fate directly on individual nuclei in a common cytoplasm, thus removing the need for transmembrane receptors, cytoplasmic signal transduction machinery and other possible events. With respect to the mode of dispersion, the prevailing model is that there is a point source of bicoid mRNA at the anterior pole of the embryo from which bicoid protein is synthesized and then diffuses posteriorly forming a gradient. Recent research suggests that a diffusion-based mechanism is sufficient to explain the exponentially decreasing gradient of bicoid protein observed (which gradually decreases as the distance from the center of production increases). However, the role of this diffusion-based mechanism on the movement of the bicoid is the subject of intense discussions. Lipshitz et al. [28] and Spirov et al. [36] suggest that the bicoid protein gradient is produced by a bicoid mRNA gradient, that is formed by active transport of a Stau-bicoid mRNA complex (i.e. bicoid mRNA in a complex with the mRNA-binding protein Staufen) through a microtubular network. Numerical solution of a the Fick’s law-based dynamical model of morphogen diffusion from the localized source of production and linear morphogen degradation shows that the morphogen gradient achieves a steady state over a broad spatial region on time periods one order of magnitude larger than the morphogen half life [4]. To demonstrate the actual participation of diffusion as a transport mechanism for the bicoid, several recent studies have been used to evaluate the bicoid’s diffusion rate [16]. The data, however, have yielded strikingly different values [15, 16]. On the initial study, injecting inert fluorescent dextran molecules into the anterior pole of the Drosophila syncytium and measuring the fluorescence intensity over 1 hrs at different spatial positions a few hundred microns from the injection Gregor et al. [15] generated time-evolution data curves which were fitted by computationally derived time courses from which the diffusion rate was calculated. Since dextran molecules of several molecular masses were used in the range of the bicoid protein molecular mass (55-57 kDa, [8]), the diffusion coefficient was uncovered using the StokesEinstein relationship, which establishes that the diffusion coefficients decrease inversely with increasing molecular radius. Thus, if the overall gradient is at a steady state and is formed by diffusion and linear degradation, based on the characteristic length of the gradient, assumed to be 100 µm, the diffusion rate was inferred to be in the order of 10 - 13 µm2 /s. However, this value does not agree with the diffusion rate inferred from the direct measurement of transgenic bicoid-enhanced green fluorescent protein levels after photobleaching at the cortical cytoplasm. The diffusion coefficient of 0.3 µm2/s, three orders of magnitude smaller than the diffusion rate doi:10.2390/biecoll-jib-2010-150 18 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de of dextran molecules and one to two orders of magnitude lower than what would be required to establish the observed bicoid gradient [16] was what was revealed. These observations are therefore inconsistent with a pure and non-generalized Fickian diffusion model, raising the open issue of how to conciliate these two disparate results. In trying to resolve this issue of disparities in data, we propose to model the bicoid protein gradient formation as a stochastic reaction-diffusion system where diffusion is a state-dependent process. By implementing this model in Redi, we attempt to point out the range of plausible bicoid protein diffusion rates. To do this, we start by assuming that the reaction-diffusion system involves the following events: • bicoid protein production; localized at the anterior pole of the embryo; • bicoid protein anisotropic diffusion; • uniform bicoid degradation. The experimental data used in this study have been recorded by T. Gregor et al. [16] in a timelapse movie (available as supplementary material of [16]) of a Drosophila Embryo expressing bicoid-GFP two-photon microscopy. The video in AVI format contains 154 frames, one frame per minute from 40 min to 194 min. The best agreement between the experiments and the simulations have been obtained for the following values of the parameters: (1) molecular mass of the bicoid protein in the range of 50-60 kDa; (2) bicoid protein production rate of order of magnitude of 10−5 min−1 ; and (3) degradation rate of bicoid protein of the order of magnitude of 10−3 min−1. The initial conditions of the spatio-temporal simulation are defined by the fluorescence values reported in the first frame of the experimental video (Figure 6 A). The reaction-diffusion system of bicoid has been simulated in a two-dimensional reaction space of 450×200 µm2 , and the cell has been fixed to 5µm, that is the radius of a nucleus [18] (1 pixel on the image measures 1µm in the physical space). Finally, the output is saved as a text file, visualized and then saved as RAW images (Figure 6 B). By using the writer HeatMapView through the command /writer:Cosbi.DiffuseSim.HeatMapViewWriter,HeatMapView a visualization of the state of the systems as a heat-colour image at each step of the simulation has been obtained (Figure 7). Due to space limitations, Figure 7 reports only 8 frames of both the experimental and the Redi simulated video that actually contain 155 frames. We can see that the Redi simulated spatio-temporal dynamics mimics quit well the observed one. The Figure 8 shows a comparison between the experimental and the simulated bicoid protein gradient. Simulations are in good agreement with the experimental observations. On the y-axis the average intensity/concentration of the bicoid, scaled in the range between 0 and 1, is reported, and on the x-axis the distance fromt he anterior pole of the embryo is reported. In order to calculate the average intensity/concentration of the bicoid corresponding to a given distance x from the tip of the bicoid on the antero-posterior axis, the image of the embryo has been divided into vertical slices 1 pixel wide, the on each slice the average pixel intensity has been calculated and normalized in the range between 0 and 1. On the graphics of Figure 8 the diffusion of bicoid away from its source can be recognized as progressive widening of the gradient doi:10.2390/biecoll-jib-2010-150 19 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de A B Figure 6: A. The initial part of the input file to Redi. The figure shows the declaration of the species (only Bicoid Protein) and the reaction (production and degradation of the protein); the parameters of the simulation are the molecular mass of Bicoid protein, the rate constants of the production and degradation reactions. The initial conditions are specified in a sort of ”matrix-like” syntax. B. The command line launched to run the Redi simulation. profile toward greater values on the x-axis. If the concentration of bicoid is mostly disposed on the anterior layer of the embryo the gradient profile is significantly different from zero for small distances x from the tip of the embryo. On the contrary, if a consistent number of bicoid molecules moved away from the anterior pole, the gradient profile will be significantly different from zero in the posterior regions of embryo (i.e. for higher values of the x coordinate). Slight discrepancies are revealed in the spatial region within 20 µm from the anterior pole of the egg, where the experimental data are most noisy and, consequently, the approximation of first derivative given in Eq. (9) is less accurate. The discrepancy between the experimental video and the simulated one has been calculated frame by frame as L1 distance. As we can see in Figure 9, it amounts to the 20-25% of the measured fluorescence. Nevertheless, the Euclidean metric to define the simulation accuracy provides an overestimate of the value of the difference between simulated and observed image. Slight effects of translation and dilatation in the simulated images are consequences of the numerical approximations and/or the propagation of errors through the steps of the algorithmic simulation procedure. To get around this, we used a distance measure based on Mahalanobis distance to take into account the translations and dilation in comparing the simulations with the experiments [6]. The Mahalanobis distance between experiments and simulations is one order of magnitude less than the Euclidean distance: 2-3% of the measured fluorescence. Figure 10 reports the average behaviour of the Mahalanobis distance for a simulation of bicoid diffusion where the molecular mass is 55 kDa. Furthermore we considered, as in [15] five nominal molecular masses in the range form 10 to 150 kDa and then launched 100 simulations for each value of the molecular mass. The best agreement, estimated by a Mahalanobis distance, between observations and simulations were obtained for simulations performed with a molecular mass between 50 and 60 kDa. The measured molecular mass of the protein is seen to drop significantly within this range (see doi:10.2390/biecoll-jib-2010-150 20 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Table 1) and indicates that the range of variability of the diffusion coefficient is between 14 and 18 µm2 s−1 , values which were able to reproduce the observed rate of the Bicoid gradient formation. Finally, we simulated the bicoid dynamics with fixed values of bicoid diffusivity (standard Fick’s law) in and out the estimated range of variability [14, 18] µm2 s−1 . We ran 100 simulations for each of the following nominal fixed values of bicoid diffusivity: 10−1 , 1, 10, 16, 30 µm2 s−1 . Table 2 shows the average Mahalanobis distance between observed and simulated dynamics. The best agreement has been achieved for a diffusivity around 10 µm2 s−1 ; nevertheless, we found that using the Fick’s law with constant diffusion coefficients does not allow to reproduce the observed movement of the bicoid protein. For instance, Figure 11 shows the simulation of the bicoid dynamics obtained using the Fick’s law with a value of the diffusion coefficient equal to 0.3 µm2 s−1 . We notice that the movement of the bicoid is significantly slow with respect to the observed one, and at the time-scale of the experiments it is not characterized by the shuttling of the bicoid in and out the nuclei of the embryo. The observed movement can be obtained with a high accuracy with our generalization of the Fick’s law. Moreover, with respect to the other diffusion-based models of the bicoid concentration profile, that are also able to reproduce the observations, our model does not need extra terms in the rate equation describing the spatio-temporal dynamics of the protein, and therefore can be considered a pure diffusive model. Namely, at the best of our knowledge, the main studies aiming at modelling the bicoid movement as a diffusive process account both of the free and the nuclear bicoid dynamics. I. Hecht et. al. and M. Coppey et al. [7, 18] reproduced the observed dynamics of the bicoid by describing with two different rate equations the spatio-temporal dynamics of free bicoid and the one of nuclear bicoid concentration. T. Gregor [15] introduced an extra negative term proportional to the local bicoid concentration in the second Fick’s law. These models reproduced the observed exponential decay of the bicoid concentration along the antero-posterior axis of the embryo. Our model and the Redi software were also able to reproduce both the temporal and spatial diffusion dynamics of the protein simulating the entire experimental movie. Bicoid protein molecular mass Average Mahalanobis distance (kDa) (order of magnitude) ≤ 10 1 40 10−1 50 to 60 10−2 70 10−1 150 1 Table 1: Order of magnitude of the discrepancy between experiments of simulations performed with five nominal molecular mass of Bicoid protein. The best agreement has been achieved for a molecular mass in the range from 50 to 60 kDa. doi:10.2390/biecoll-jib-2010-150 21 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Figure 7: Experimental video frames (black figures on the left) and Redi simulated video frames (red heat-color maps on the right). The simulated images are an average of 100 stochastic simulations. doi:10.2390/biecoll-jib-2010-150 22 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Figure 8: Average Bicoid protein fluorescence profiles along the antero-posterior axis at different instants of time. The black curve is the experimental profile and the red is the curve obtained from the Redi simulation. Simulations are in good agreement with the experimental observations. Discrepancies are visible mostly as a slight shift of the simulated curve with respect to the measured one in the first minutes of the simulation and within a distance of 20 µm from the anterior pole of the egg. In this region the experimental video is noisy and, consequently, the first approximation of the derivative of the concentration in Eq. (9) is not accurate enough. 6 Discussion The anisotropic stochastic diffusion model presented in this paper reproduces the spatio-temporal dynamics of the bicoid observed with direct photobleaching measurements by T. Gregor et al [16]. The model is a generalization of the Fick’s law in which the diffusivity is function of the time and on the local concentrations. The initial conditions are defined by values of intensity/doi:10.2390/biecoll-jib-2010-150 23 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Figure 9: Euclidean distance between experimental frames and Redi simulated frame. concentration on the first snapshot of the time-lapse movie of the Drosophila embryo recorded by T. Gregor [16]. The bicoid gradient is one of the best studied morphogen gradients. During the oogenesis its mRNA is deposited at the anterior pole of the egg and is translated soon after fertilization. Antibody staining in fixed embryos demonstrates that bicoid concentration decays exponentially with distance from the anterior pole [8, 15, 16]. The observed exponential shape is consistent with the gradient being formed by a balance of localized synthesis, diffusion, and spatially uniform degradation, which are normally referred to as the SSD model. Important parameters in an SSD model are the diffusion coefficient and the bicoid lifetime in the cytoplasm. doi:10.2390/biecoll-jib-2010-150 24 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Figure 10: Time behavior of average Mahalanobis distance between experimental and simulated spatio-temporal dynamics of Bicoid protein gradient. Bicoid diffusion constant (µm2 s −1 ) 10−1 1 10 16 30 Average Mahalanobis distance (order of magnitude) 1 10−1 10−2 10−2 102 Table 2: Order of magnitude of the discrepancy between experiments of simulations performed with five nominal constant values of bicoid diffusivity. The best agreement has been achieved for a diffusivity around 10 µm2 s −1 . Figure 11: Bicoid spatio-temporal dynamics simulated with Fick’s law with constant diffusion coefficient equal to 0.3 µ m 2 s−1 . The kinetics results to be much slower than the observed one and it no shuttling movement of bicoid in and out the nuclei are reproduced in the time interval from 10 to 150 min. Various variants of SDD model have been proposed to fit the experimental data and estimate the kinetic parameters of the bicoid dynamics [7, 15, 16, 18], as a simple diffusion/degradadoi:10.2390/biecoll-jib-2010-150 25 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de tion model turned out to be inconsistent with the experimental observations. A Fick’s diffusion √ and a constant degradation rate δ and a diffusivity D would lead to a length scale λ = DT where T is the time interval during which bicoid has diffused and λ is the length over which the profile decays from its maximal value (in the anterior pole of the egg) to zero. However this calculation lead to length scales that are much smaller than the observed ones. Furthermore, simple SDD models are unable to capture the rapid nucleocytoplasmic shuttling observed in the live-imaging experiments of T. Gregor et al.. M. Coppey et al. [7] and I. Hecht [18] proposed to similar models of diffusions and nuclear trapping able to estimate the observed length scale, steady profile and scaling between embryos of different lengths. Here we briefly review them and compare them with our approach. Consider a one-dimensional model of the embryo of length L. Bicoid molecules are produced at a constant rate kprod at the anterior and of the embryo (x = 0); the posterior of the embryo (x = L) is impermeable to bicoid. For the dynamics model, a 1D model of teh embryo cortex with no boundary conditions at either end is used. This 1D model is a projection of the embryo onto the antero-posterior axis. This axis is discretized using a discretization parameter ∆x and nuclei are positioned at xn . In [7] the embryo is modeled as a homogeneous medium, where nuclei are uniformly distributed. Bicoid can exist in two states: ”out” when it is out of the nucleus, and ”in” when it is into the nucleus. The transition between ”in” and ”out” are modeled by first order reactions [7, 18] (Figure 12). Based on this, the dynamics of the bicoid concentration is described by the system of equations (41)-(44). ∂Bin = kin Bout − kout Bin − δin Bin ∂t ∂ 2 Bout kin ∂Bout =D − δout Bout − Bout δ(x − xn ) + kout Bin 2 ∂t ∂x A∆x ∂Bout D = −kprod ∂x ∂Bout =0 D ∂x (41) (42) (43) (44) The concentration of bicoid in the nucleus is determined by the flux of extranuclear free bicoid kin into the nucleus, by the flux of intranuclear bicoid kout coming out of the nucleus, and by the nuclear degradation rate δin . The extranuclear concentration of bicoid is determined by the diffusion motion (described by the second Fick’s law in Eq. (42)), by the extranuclear degradation rate δout , by the flux of bicoid kin going into the nucleus. In this 1D model, A∆x is a geometric factor arising from the projection onto the antero-posterior axis of the embryo; A is the area of the cross section of the embryio that in this model is supposed to be constant throughout the embryo [18]. The study of T. Gregor et al. [16] showed that there is a dynamics equilibrium between transport into the nucleus and a loss, either via outward transport or intranuclear protein degradation. This leads to Bin =K Bout where K is the equilibrium constant for the nucleocytoplasmic shuttling. If the nuclei are uniformly distributed in the embryo, as assumed in [7], K is a function of time, but does not depend on the spatial coordinates. From this, the local concentration of nuclear and extranuclear bicoid can be expressed as the total concentration of bicoid Btot (x, t) = Bin (x, t) + Bout (x, t), doi:10.2390/biecoll-jib-2010-150 26 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de Figure 12: Bicoid exists in two states: freely diffusing Bout , and immobile/nuclear Bin . The transitions between the two states are described by first order processes [7, 18]. The forward nuclear trapping rate constant is ki n, and the reverse process rate constant is kin . so that Btot K Bin = Btot . 1+K 1+K Consequently, adding Eq. (41) to Eq. (42) leads to Bout =  
1 ∂Btot 1 D ∂ 2 Btot δ(x − xn ) − = k − 1 B δ − + δ K Btot (45) in tot out in 2 ∂t A∆x 1+K {z } |1 + K{z ∂x } |1 + K | {z } diffusion nuclear trapping degradation The Eq. (45) consists of three terms: the first models the diffusion, the second models the nuclear trapping, and finally the third models the degradation. Hecht’s et al. model [18] includes also a space-dependent production function P (x), and a cytoplasmic flow defined as ~νflow (x) · ∇B. This term describes the coupling pf the diffusion to the fluid flow where ~νflow (x) is the fluid velocity. The existence of the cytoplasmic flow could explain the motion of the nuclei in the viscous cytoplasm. Our approach to bicoid dynamics includes stochastic production, degradation, and diffusion of the local concentration. All these processes are modelled as first order processes. The diffusion coefficient is a function of the spatial coordinates, so that, based on the local conditions, it could change box-to-box onto the spatial grid partitioning the physical space. The size of the boxes of the spatial grid given by Eq. (32), is equal to the average radius of the nuclei, as a consequence a box can contain only one nucleus. The spatial distribution of nuclei is not uniform, therefore, some boxes contain a nucleus, while some others do not. When the diffusion of bicoid diffusion take place between to adjacent boxes both containing a nucleus, molecules of bicoid are exchanged between the neighboring nuclei. Otherwise, when the diffusion takes place between to adjacent boxes of which only one contains a nucleus, molecules of bicoid are exchanged between the nucleus and the extranuclear space. Since the diffusion coefficient is space-dependent, if the local conditions (concentration, friction coefficient, chemical activity) in the two boxes are different, the diffusion coefficient will take a different values in the two boxes. The Figure 13 shows a situation in which a box is occupied by a nucleus and the next one is not. The diffusions of the bicoid in and out of the nucleus are modelled as first order events having rate constant kin = Din /l2 and kout = Dout /l2 , respectively, where     kB T Bin ∂γin kB T Bout ∂γout Din = 1+ , Dout = 1+ fin γin ∂Bin fout γout ∂Bout doi:10.2390/biecoll-jib-2010-150 27 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de where γin (γout ), fin (fout ) are the chemical activity and the frictional coefficient defined by the local condition in and out of the nucleus, respectively. Figure 13: The diffusion event from box to box of the grid representing the physical space is modelled as a first order reaction whose rate constants is proportional to the diffusion coefficient. In our model the diffusion coefficient is space-dependent, so that it has different values in different boxes depending on the local values of bicoid concentration, frictional coefficient and chemical activity. Therefore, our approach based on a space-dependent model of diffusion coefficient and on a first order representation of the diffusional movement of the bicoid molecules, is capable to simulate both the exchange of molecules between neighboring nuclei as well as the passage of molecules in and out of the nuclei. Moreover, our simulation framework is stochastic, i. e. degradation, production and diffusion are stochastic events selected on the basis of their propensities [14]. 7 Conclusions and future directions The motivation for this study was inspired by recent wet-lab experiments confirming that constant diffusion coefficients are rather more the exception than the rule in living cells and, more generally in biological tissues. We have presented in this report a model that is able to represent the diffusion of non-charged molecules, in which the diffusion coefficients are not constant with respect to the time and space. To show this, we described a model of a state-dependent diffusion able to simulate the physical reality of biochemical reaction-diffusion processes. We designed and built the software tool Redi to implement our model of diffusion in the framework of stochastic simulation of reaction-diffusion systems. Then we presented the results of the Redi simulations on the case study of the spatio-temporal dynamics of the bicoid protein. Using Redi we were able to show that state-dependent diffusion plays a role in the bicoid transport and were able to pinpoint a molecular mass for bicoid that reproduces the measured bicoid movement in the Drosophila embryo with the highest accuracy. Regarding this cases study, new models and simulations are currently under development to determine if the stochastic reaction diffusion system is able to simulate the spatio-temporal dynamics of the bicoid protein with a non-localized source. Recent studies indicate that the bicoid mRNA itself if distributed in a gradient and not localized at he pole [36]. Unlike previous works [5, 19, 34], this model provides a theoretical derivation of the molecular origins of the parameters, determining the time-behaviour of the diffusive phenomena in non-homogeneous media. The model of diffusion proposed in this report is general and was doi:10.2390/biecoll-jib-2010-150 28 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de successfully applied to simulate the spatial dynamics of any molecule in non-homogeneous media. Some case-studies where spatial effects are relevant which we have also successfully simulated with Redi are: tubulin diffusion in cytoplasm [27], and chaperone-assisted protein folding [26]. Future works will consist in further refinement of the procedure to take into account the chemistry and physics of biological transport phenomena. Some future directions will consist of a more accurate calculation of the second virial coefficient for biomolecules, especially for proteins. The use of the Lennard-Jones potential is a good approximation of the molecular interaction, but is a drawback in describing protein-protein interactions of water molecules must be included explicitly [32], complicating the computational task. The condition of solvated molecules is reflected also in the expression of the concentration-dependence of frictional coefficient, that will need to be modified accordingly. More generally the cellular environment is a crowded solution. Namely, the cellular environments are packed with other biomolecules and this crowdedness may affect the stability and aggregation rates of proteins inside cells [21, 22, 13, 29, 30]. Unlike in typical biochemical experiments in which the proteins of interest are purified and diluted, the living cell is crowded with a wide variety of other proteins and macromolecules which generally occupy 20-30% of the total cell volume. This percentage is called excluded volume. The effects imposed by the excluded volume are caused by the volume excluded by the “inert” macromolecules, are called macromolecular crowding effects and those macromolecules are called crowding agents. We are currently extending the present study to develop a model whose simulations are able to support the investigation of excluded volume effects on the protein diffusion and folding. Other improvements are currently under development, and rely on the stochastic algorithm of Redi. It can be incorporated with the time extension of Gillespie algorithm, that the authors developed for the simulation of process-based models, to treat rate coefficients depending on time [24, 25]. References [1] P. S. Agutter and D.N. Wheatley. Random walks and cell size. BioEssays, 22:1018–1023, 2000. [2] P.S. Agutter, P.C. Malone, and D.N. Wheatley. Intracellular transport mechanisms: a critique of diffusion theory. J. Theor. Biol., 176:261–272, 1995. [3] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter. Molecular biology of the cell. Garland Science, 4th ed. edition, 2003. [4] S. Bergmann. Pre-steady state deconding of the bicoid morphogen gradient. PloS Biol., 5:232, 2007. [5] D. Bernstein. Simulating mesoscopic reaction-diffusion systems using the gillespie algorithm. Physical Review E, 71, April 2005. [6] X. Chen and T. J. Cham. Discriminative distances measures for image matching. In Procs of Int. Conference on Pattern Recognition, Cambridge, England, 3:691–695, 2004. doi:10.2390/biecoll-jib-2010-150 29 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de [7] M. Coppey, A. M. Berezhkovskii, Y. Kim, A. N. Boettiger, and S. Y. Shvartsman. Modeling the bicoid gradient: diffusion and reversible nuclear trapping of a stable protein. Developmental Biology, 312:623–630, 2007. [8] W. Driever and C. Nuessein-Volhard. A gradient of bicoid protein in drosophila embryos. Cell, 54(1):83–93, 1988. [9] J. Elf, A. Doncic, and M. Ehrenberg. Mesoscopic reaction-diffusion in intracellular signaling. Fluctuation and noise in biological, biophysical and biomedical systems. Procs. of SPIE, 5110, 2003. [10] J. Elf and M. Ehrenberg. Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases. Syst. Biol., 1(2), December 2004. [11] M. B. Elowitz, P. E. Wolf M. G. Surette, J. B. Stock, and S. Leibler. Protein mobility in the cytoplasm of escherichia coli. J. Bacteriol., 181:197–203, 1999. [12] D. Fusco, N. Accornero, B. Lavoie, S. Shenoy, J. Blanchard, R. Singer, and E. Bertrand. Single mrna molecules demonstrate probabilistic movement in living mammallian cells. Curr. Biol., 13:161–167, 2003. [13] G. Yang G. Ping and J. M. Yuan. Depletion force from macromolecular crowding enhances mechanicsl stability of protein molecules. Polymer, 27:2564:2570, 2006. [14] D.T. Gillespie. Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry, 81(25):2340–2361, December 1977. [15] T. Gregor, W. Bialek, R. R. de Ruyter van Stevenick, D. W. Tank, and E. F. Wieshaus. Diffusion and scaling during early embryonic pattern formation. PNAS, 102(51):18403– 18407, 2005. [16] T. Gregor, E. F. Wieshaus, A. P. McGregor, W. Bialek, and D. W. Tank. Stability and nuclear dynamics of the bicoid morphogen gradient. Cell, 130:141–152, 2007. [17] S.E. Harding and P. Johnson. The concentration dependence of macromolecular parameters. Biochemical Journal, 231:543–547, 1985. [18] I. Hecht, W. J. Rappel, and H. Levine. Determining the scale of the bicoid morphogen gradient. PNAS, 106(6):1710–1715, 2009. [19] S. A. Isaacson and C. S. Peskin. Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations. SIAM Journal of Scientific computing, pages 47–74, 2006. [20] E. R. Kandel. The molecular biology of memory storage: a dialogue between genes and synapses. Science, 294:1030–1038, 2001. [21] A. R. Kinjo and S. Takada. Effects of macromolecular crowding on protein folding and aggregation studied bu density functional theory: statics. Physical Review E, 66:031911: 1–9, 2002. doi:10.2390/biecoll-jib-2010-150 30 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de [22] A. R. Kinjo and S. Takada. Effects of macromolecular crowding on protein folding and aggregation studied by density functional theory: Dynamics. Physical review. E, 66(5):051902.1–051902.10, 2002. [23] K.J. Laidler, J.H. Meiser, and B.C. Sanctuary. Physical Chemistry. Houghton Mifflin Company, 2003. [24] P. Lecca. A time-dependent extension of gillespie algorithm for biochemical stochastic π-calculus. Proceedings of the 2006 ACM symposium on Applied computing, 2001. [25] P. Lecca. Simulating the cellular passive transport of glucose using a time-dependent extension of gillespie algorithm for stochastic π-calculus. International Journal of Data Mining and Bioinformatics, 1(4):315–336, 2007. [26] P. Lecca and L. Dematté. Stochastic simulation of reaction-diffusion systems. Int. F. og Medical and Biological Engineering, 1(4):211–231, 2008. [27] P. Lecca, L. Dematté nd M. Lecca, and C. Priami. Stochastic modelling of diffusion systems. video image simulation of tubulin diffusion in cytoplasm: a case study. In II Eccomas thematic conference on computational vision and medical image processing, 2009. [28] H. D. Lipshitz. Follow the mrna: a new model for bicoid gradient formation. Nature reviews, 10:509–512, 2009. [29] A. P. Minton. Molecular crowding: analysis of effects of high concetrations of inert cosolutes on biochemical equilibria and rates in terms of volume exclusion. Methods Enzymol., 295:127:149, 1998. [30] A. P. Minton. The influence of macromolecular crowding and macromolecular confinement on biochemical reactions in physiological media. J. Biol. Chem., 276:10577:10580, 2001. [31] M. J. Moran and H. N. Shapiro. Fundamentals of Engineering Thermodynamics. Wiley, 3th ed. edition. [32] B. L. Neal, D. Asthagiri, and A. M. Lenhoff. Molecular origins of osmotic second virial coefficients of proteins. Biophysical Journal, 75, 1998. [33] C. Priami. Algorithmic systems biology. Communications of the ACM, 52(5):80–88, 2009. [34] C. J. Roussel and M. R. Roussel. Reaction-diffusion models of development with statedependent chemical diffusion coefficients. Progress in Biophysics & Molecular Biology, 2004. [35] A. Solovyova, P. Schuck, L. Costenaro, and C. Ebel. Non ideality of sedimantation velocity of halophilic malate dehydrogenase in complex solvent. Biophysical Journal, 81:1868–1880, 2001. [36] A. Spirov, K. Fahmy, M. Schneider, E. Frei, M. Noll, and S. Baumgartner. Formation of the bicoid morphogen gradient: an mrna dictates the protein gradient. Development, 136:605–614, 2009. doi:10.2390/biecoll-jib-2010-150 31 Journal of Integrative Bioinformatics, 7(1):150, 2010 http://journal.imbio.de [37] M. P. Tombs and A. R. Peacocke. The Osmotic Pressure of Biological Macromolecules. Monograph on Physical Biochemistry, Oxford University Press, 1975. doi:10.2390/biecoll-jib-2010-150 32