Modeling Reaction Kinetics in Low-dimensional Environments with Conformon P Systems: Comparison with Cellular Automata and New Rate Laws

Modeling Reaction Kinetics in Low-dimensional Environments with Conformon P Systems: Comparison with Cellular Automata and New Rate Laws Pierluigi Frisco1, Ramon Grima2 1School of Mathematics and Computer Sciences Heriot-Watt University, Edinburgh, UK [email protected] 2School of Biological Sciences University of Edinburgh, UK [email protected] Summary. Recently it has been shown that simulations of complex biological systems using conformon P systems and cellular automata do not necessarily give the same predictions. To further elucidate these differences we simulate a simple model of intracellular reactions involving a single bimolecular reaction occurring on a biological membrane using conformon P systems. We find that the predictions broadly agree with results from both the theory of random walks in low-dimensional environments and with previously published simulations using cellular automata. Moreover, a re-analysis of the data enables us to deduce novel rate laws for the kinetics of reactions occurring on biological membranes. 1 Introduction A recent publication [3] reported that simulations of HIV dynamics differ in their results according to the simulation platform used. In particular it is found that cellular automata (CA) models produce qualitatively correct dynamics only for a narrow range of rule probabilities and for particular initial conditions whereas conformon P (cP) models [2] derived from the CA model display significant robustness of qualitatively correct dynamics over a wide range of conditions. Presently the reasons for these differences are not understood. The complexity of the system under study precludes a rigorous analysis of these discrepancies. In this paper we consider a much simpler biological process at the base of an simpler model. For such model its rigorous analytical results are known for some cases. The paper is divided as follows: in Section 2 we describe the biology behind the model and its implementation using CA and cP, in Section 3 we analyze the 134 P. Frisco, R. Grima data generated by cP, compare with theory and CA results and also deduce some new biologically relevant kinetic laws and in Section 4 we draw some conclusion. 2 A Simple Biological Model: Implementation Using CA and cP There is growing evidence of the importance of reaction kinetics for the structural organization of the intracellular environment, which is far from the homogeneous, well mixed solution typical of in vitro experiments. A high degree of molecular crowding as well as the presence of indigenous obstacles in cellular media have important consequences in the physico-chemistry of the cell. The consequences of this in some process are only now becoming to be more generally understood. One of these consequences is that it is not clear what are the rate laws governing reactions occurring in vivo [7]. To tackle this problem biochemists have been using various computational frameworks to extract rate laws or empirical reaction equations from direct numerical simulations. Among these approaches, simulations based on CA are the most popular (see, for example, [5]). The biological process we considered in our investigations regards biochemical reactions occurring on cell membranes. It in known that about half of the proteins inside cells are membrane-associated [8] and thus biochemical reactions must necessarily function within the constraints imposed by the two-dimensional environment of the biological membrane. Prominent examples of such reactions are those involving enzymes called lipases which play key roles in fat metabolism and digestion and which occur on two-dimensional interfaces rather than in threedimensional solution. The simplest model of such dimensionally-restricted reaction kinetics consists of two types of particles, denoted with Aand B, which perform random walks on a two-dimensional plane and which upon encounter react with some probability and produce a single new inert particle C. This mathematical construct represents the physical process of the reaction of two molecules of two different types which normally perform Brownian motion (modelled by the random walks) and which react upon encounter to form some new product molecule [4]. Such elementary reactions form the backbone of all biochemical reaction networks, independent of their complexity and are particularly ideal for a comparison between CA and cP models because of the existence of rigorous analytical results from the theory of random walks in low-dimensions. The biological process indicated above simplifies the biological membrane to a homogeneous quasi-two dimensional environment. In reality it is found that the heterogeneous micro-structure of the membrane significantly hinders the free diffusion of molecules on its surface. In particular it is known that transmembrane proteins (denoted with Bin the above) impose relatively static barriers to the smaller and more mobile molecules (denoted with Ain the above). This is due because transmembrane proteins are anchored to the cytoskeleton of the cell. These obstacles are incorporated in the models considered by us by making some parts of the plane inaccessible to particle motion. Modeling Reaction Kinetics with Conformon P Systems 135 A general CA model which describes both cases above (with and without obstacles to particle motion) has been described in [4]. The algorithm is the follows. Initially, particles of two different kinds Aand Bare uniformly distributed on a two-dimensional lattice with unit spacing and periodic boundary conditions (a torus). Particles Acan move, while particles Bare static. One Aand one Bparticle can react when in the same location of the lattice and produce one particle C. Particles Care static and inert. Some of the tests considered another type of static and inert particle, an obstacle, uniformly distributed in the lattice in the initial configuration. At each time step, a particle of type Aor Bis randomly chosen and either moved or subject to a reaction according to the following: the particle can move from the location in the lattice in which it is to a randomly chosen neighbor location only if the chosen neighbor location does not contain any other particle of the same kind or an obstacle; if instead the chosen neighbor lattice location contains a particle of the other kind (that is, if Ais subject to be moved, then Bis the other kind; if Bis subject to be moved, then Ais the other kind), then the particle can react with the particle of the other kind with probability P. If this occurs both particles A and Bare removed and a Cparticle is placed in the chosen neighbor location, otherwise nothing occurs. It is important to note that the algorithm does not allow more than one particle of any type to be in the same location of the lattice, thus enforcing a hard-sphere molecular repulsion. The above two steps are repeated ntot(t) times, where ntot(t) is the current number of distinct particles on the lattice (excluding obstacles) at time t. After one such sequence the time is incremented by one. The simulations are performed with two different lattice types: square (von Neumann) and triangular neighborhoods. Models of cP systems have been derived by the just described CA model. Particles of type Aand Bhave been modeled with [A, 1] and [B, 1] conformons, respectively. Their eventual interaction (with probability P) creates [B, 2] conformons representing the Cparticles. A lattice location of the CA has been modeled with a (membrane) compartment in the cP model. The presence of obstacles has been modeled with compartments in the lattice having no incoming edge (in this way no conformon could move in these compartments). The simulations have been performed using the cP simulator available from [9] modifying it in a way that interaction rules have priority on passage rules. Moreover, ad hoc programs to create the lattices and to analyze the data have been also used. These programs can be requested to the authors. The cP models and the simulations are such that more than one particle can be at the same time in one compartment. 3 Data Analysis Data produced by the simulations consisted of the number of Aand Bparticles as a function of time. For each set of parameter values, ten independent simulations 136 P. Frisco, R. Grima were run on 300×300 lattices until the time (denoted by T) at which there remained only 1% of the minority particle species (that is, the one whose initial concentration is the smallest, which in our case are particles of type A[4]). The majority species is static, the minority species is mobile and their probability of interaction is P. We run the tests with two different values of P: 1 and 0.1. We consider two cases: presence and absence of obstacles uniformly distributed throughout the lattice so to occupy at most 0.4 of the available lattice locations (denoted by [O]=0.4). The data from these ten copies were then averaged to reduce the inherent noise, yielding an array of values [A(t)],[B(t)], t, t = 0...T. The simulation data for the absence of obstacles case can be directly compared to rigorous theoretical results from the theory of random walks [1]. Here it is shown that in a two-dimensional space where the majority particle species is immobile (B) and minority species is mobile (A) and the reaction probability is 1, then for sufficiently long time runs, the rate of change of the concentrations of mobile particles is described by the effective ordinary differential equation: d[A(t)]/dt ∼ −t−1/2[A(t)][B(t)],(1) where [A(t)] denotes the total number of particles of type Aat time tdivided by the total number of lattice locations defining the two-dimensional space. Thus the most basic test of our cP simulations is to use the data obtained for the case of no obstacles and P= 1 to extract the time exponent in the above equation. The method used to obtain this exponent is the one reported in [4] where it is shown that for general differential equations of the type d[A(t)]/dt ∼ −t−(1−p)[A(t)][B(t)], the exponent pis equal to the gradient of the graph of G=Log[−Log(B0[A(t)]/(A0(B0−A0+ [A(t)])))] versus Log(t). Figure 1, bottom curve, shows the variation of the slope (that is, p) with time for the just indicated cP simulations. It is found that p= 0.6. This implies that the ordinary differential equation satisfied by the cP simulation data is d[A(t)]/dt ∼ −t−0.4[A(t)][B(t)]. This result is fairly close to the rigorous theoretical value given above and also agrees with previous CA simulations giving p= 0.5. Figure 1, above curve, shows the results of the simulations for P= 0.1. No rigorous theory exists for this case, but CA simulations [4] give p∼0.85 whereas with cP simulations we obtain p∼0.92. Hence for the case of no obstacles, for both high and low values of reaction probability P, the results of CA and cP are in good quantitative agreement, though there is a consistent tendency of the exponent for cP to be slightly larger than that of CA. The latter discrepancy could be due to the fact that CA simulations impose the condition that only one particle is allowed at a site whereas cP simulations make no such assumption. The lack of such an assumption would necessarily imply a larger amount of “particle mixing” inside each spatial element of cP simulations which from physical considerations [4] would necessitate the exponent pto be closer to one, as observed. To test this hypothesis we developed a cP model in which at most one particle per compartment is allowed. Figure 2 shows the curve comparing the data obtained by the cP models in which more than one particle and at most one particle per Modeling Reaction Kinetics with Conformon P Systems 137 Fig. 1. Variation of the slope with no obstacles in the cP model compartment is allowed. No discernible differences are observed between the two models. This implies that either the effect occurs only in lower dimensions or the possibility to have more than one particle in the same compartment is not the reason for the small discrepancies between CA and cP simulations. There are no rigorous analytical results for the case in which obstacles are present. Anyhow, it has been traditionally assumed that the dynamics in this case would be captured by an effective ordinary differential equation as (1), but with a time exponent pwhich varies somewhere between 0 and 1. This is often referred to as fractal kinetics [5, 7]. However, in [4] it is shown that this is not the case. It is found from CA simulations that the slope is not constant but varies considerably with time and apparently does not approach a constant value in the limit of long time runs. Our cP simulations also confirm this result (Figure 3), once again showing no evidence of a discrepancy between CA and cP simulations. In the present paper we go one step beyond the work reported in [4] and, for the case in which obstacles are present, we find a new effective ordinary equation which captures the dynamics of the reaction. It can be shown [4] that the solution of an ordinary differential equation of the type d[A(t)]/dt ∼ −k(t)[A(t)][B(t)] for long time runs is of the form: [A(t)] ∝exp·−(B0−A0)Zt 0 k(s)ds¸.(2) 138 P. Frisco, R. Grima Fig. 2. Curves comparing results from cP simulations. Blue and pink curves for Aand Bwhen at most 1 particle per compartment is allowed. Yellow and cyan curves for Aand Bwhen more than one particle per compartment is allowed. The graph shows that the variation of particle concentration with time is independent of the one-particle constraint. It is also known that k(t) = dS/dt where Sis the mean number of distinct lattice locations visited by a particle moving in a random walk [5]. It is found that S∼ tds/2for long time runs in a fractal space of spectral dimension ds. This would imply (for long time runs) the dynamics to follow an effective equation of the form d[A(t)]/dt ∼ −t−(1−ds/2)[A(t)][B(t)] which implies a constant time exponent p. However note that to arrive at this conclusion one implicitly assumes that the long time regime is being observed. Actuality one may only observe the early and intermediate time regimes since the simulation halts after 99% of the particle A has been consumed. Inspired by the theoretical results reported in [4], we surmise that the intermediate time scaling for Swould be of the general form: S∼tαLn(1/Ln(tα)) where the exponent αis introduced to take into account the heterogeneity of space imposed by the presence of obstacles. Interestingly, it is found that the cP data is in good agreement with this conjectured law, see Figure 4. Thus our simulations and data analysis suggest a new kinetic equation for describing bimolecular reactions in obstacle-ridden low dimensional media, namely d[A(t)]/dt ∼ −k(t)[A(t)][B(t)] with k(t) = ∂/∂ttαLn(1/Ln(tα)) instead of the customary k(t) = tα. Modeling Reaction Kinetics with Conformon P Systems 139 Fig. 3. Gvalue for the cP tests with more than one particle per compartment and no obstacles 4 Conclusion In this paper we report a preliminary study aiming to understand what kind of biological processes are better fit to be modeled with CA or with cP. In particular, we focused on the possibility offered by cP to model the presence of more than one particle in a compartment. From the tests we run we conclude that this possibility does not always make a difference in the obtained results. Some differences between the results obtained by similar CA and cP models occurs only if obstacles, locations in the lattice limiting the mobility of the particles, are present. Anyhow, the found differences are not yet sufficient to draw general conclusions. One line of further research is to compare the fluctuations from the average values in CA and cP simulations. Acknowledgements: R. Grima gratefully acknowledges support from SULSA (Scottish Universities Life Sciences Alliance). 140 P. Frisco, R. Grima Fig. 4. Results for the cP tests with more than one particle per compartment and obstacles References 1. M. Bramson, J.L. Lebowitz: Asymptotic behavior of densities in diffusion-dominated annihilation reactions. Phys. Rev. Lett., 61 (1988), 2397. 2. P. Frisco: The conformon-P system: A molecular and cell biology-inspired computability model. Theoretical Computer Science, 312 (2004), 295. 3. P. Frisco, D.W. Corne: Dynamics of HIV infection studied with cellular automata and conformon-P systems. Biosystems, 91 (2008), 531. 4. R. Grima, S. Schnell: A systematic investigation of the rate laws valid in intracellular environments. Biophys. Chem., 124 (2006), 1. 5. R. Kopelman: Rate-processes on fractals: theory, simulations and experiments. J. Stat. Phys., 42 (1986), 185. 6. H. Larralde, P. Trunfio, S. Havlin, H.E. Stanley, G.H. Weiss: Number of distinct sites visited by N random walkers. Phys. Rev. A, 45 (1992), 7128 . 7. S. Schnell, T.E. Turner: Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws. Prog. Biophys. Mol. Biol., 85 (2004), 235. 8. P. Reis, K. Holmberg, H. Watzke, M.E. Leser, R. Miller: Lipases at interfaces: A review. Adv. Coll. Int. Sci., 2008, doi:10.1016/j.cis.2008.06.001. 9. The P Systems Webpage: http://ppage.psystems.eu/.