Metabolic Algorithm with Time-varying Reaction Maps

Metabolic Algorithm with Time-varying Reaction Maps Luca Bianco, Federico Fontana, Vincenzo Manca University of Verona Department of Computer Science strada Le Grazie, 15 37134 Verona, Italy E-mail: {bianco,fontana}@sci.univr.it,[email protected] Summary. A symbolic-based approach to modelling biochemical processes and cellular dynamics is likely to turn useful in computational biology, where attempts to represent the cell as a huge, complex dynamic system must trade with the linguistic nature of the DNA and the individual behavior of the organelles living within. The early version of the metabolic algorithm gave a first answer to the problem of representing oscillatory biological phenomena, so far being treated with traditional (differential) mathematical tools, in terms of rewriting systems. We are now working on a further version of this algorithm, in which the rule application is tuned by reaction maps depending on the specific phenomenon under consideration. Successful simulations of the Brusselator, the Lotka-Volterra population dynamics and the PKC activation foster potential applications of the algorithm in systems biology. 1 Introduction Symbolic rewriting has traditionally been used to study and classify formal languages [18]. It was some years after Chomsky’s fundamental discoveries that rewriting systems began to be applied to the study of the growth of some simple organisms and to the analysis of biological structures [9, 15]. These early applications of rewriting to practical case studies taken from the real world demonstrated the potential ability of a properly defined formal construct to represent, in principle, the development of at least some biological species. Such constructs, in fact, move step by step toward the definition of a language/structure until their computation terminates, hence their application to species in development emphasized the possibility for formal systems to figure out not only classes of languages, but also the paths along which their final structure takes form during the system evolution. Recently, a research line has started which focuses on the rewriting system dynamic activity instead of its expressive power evaluated in terms of language 44 L. Bianco, F. Fontana, V. Manca types [20, 3, 11]. This line has been stimulated in an attempt to capture, by means of rewriting systems, the dynamics of a biochemical process. In this attempt a novel construct known as P system has come useful, provided its capability to represent several structural aspects of the cell along with many intraand extra-cellular communication mechanisms [16, 4]. Such a dynamic perspective on rewriting employing P systems has already led to alternative representations of different biological dynamics [1, 19] and to new models of important pathological processes [3, 14]. In [3, 11] we have started to develop a metabolic algorithm that introduced a new perspective in the rewriting mechanism of P systems: 1. rules are not applied to objects. Rather, they are applied to populations of objects; 2. rules are specified along with reactivities. Every reactivity denotes the ability of the corresponding rule to compete against other rules in capturing part of a population, on which the reaction is performed. We have gone further in this perspective, by associating every reactivity to a map that depends on the state of the system. Moreover we have added a strategy for partitioning the objects in the system at every transition, depending on the relative magnitude of every reactivity. The performance shown by the metabolic algorithm in the simulation of wellknown biochemical models, such as the Lotka-Volterra population dynamics [10, 21], the BZ chemical reaction [5], and the PKC activation process [2], fosters potential applications in critical open problems dealt with by systems biology [7]. Simulation in progress are confirming the effectiveness of the algorithm in modelling even more complex biochemical processes, such as those that evoke circadian rhythms in living bodies [8]. 2 Metabolic Algorithm As we have told in the introduction, the metabolic algorithm is built on P systems. For the sake of simplicity here, and in the following, we hypothesize that this system is made of just one membrane. In Section 5 we will briefly discuss the formal extensions needed to cope with more membranes. To provide our algorithm with flexibility we will guarantee the fulfillment of the two following principles at any transition of a P system Π, working on the alphabet A={X, Y, . . . , Z}and provided with rules r, s, . . . , w ∈R: •increasing the activity of a rule implies a proportional decrease of the rewriting activity of other rules sharing the same symbols. This condition reflects the concurrency among rules over a finite set of objects in the system; •the applicability of rules is limited by those objects, whose availability in the system is low. This condition reflects a constraint on resources. Metabolic Algorithm with Time-varying Reaction Maps 45 In the following we will reformulate these two points in quantitative terms. In the early version of the metabolic algorithm we had postulated that proper reactivity constants affected the rewriting activity of every rule, respectively, in a way that a larger reactivity constant defined a higher rewriting activity of the corresponding rule. Then, this activity was properly limited by defining a population on which a rule r:αr→βr, transforming the string αr∈ A∗into a new string βr∈ A∗, could be applied during a transition depending on the number of objects available in the system immediately before that transition. The new version of the metabolic algorithm requires, firstly, to recognize the state of the system. This state is used inside so-called reaction maps which generalize the former reactivity constants into time-varying functions. Once we have such maps at hand we will let the rules work according to the relative reactivity expressed by every map, meanwhile limiting this power to avoid over-consumption of the objects in the system. Finally, a simple stochastic method will be adopted to decide how to treat individual objects in the system, for which the procedure described so far does not take a definite decision. More in general, introducing stochastic properties in an algorithm can turn out to be particularly desirable when the dynamics is highly influenced by few molecules [6]. 2.1 State of the System In classical dynamic systems the values assumed along time by every variable usually form the state of the system. In a similar way, here we postulate that at every discrete time tthe number of objects of each type is well defined for every membrane. Formally, the state of the system at time tis identified by a function qt:A −→ N,(1) where Ais the alphabet of the P system. For instance, qt(X) gives the amount of objects Xavailable in the system at time t. Note that we will usually omit to denote t, except for those cases in which specifying the time turns out to be convenient. The state, hence, can be read by applying qtto every symbol of the alphabet. The set of all states assumed along time by the system is expressed by: Q={qt|t∈N}.(2) This set, then, contains the complete information on the evolution of the system. Further insight on the state is not possible since qtis the only probe we can use to observe it. 2.2 Reaction Maps As opposite to the early metabolic algorithm [3] in which the reactivity constants had a direct (and time invariant) role on rewriting, here we generalize such constants into reaction maps, one for each rule, in a way that every reaction map gives 46 L. Bianco, F. Fontana, V. Manca the reactivity that the corresponding rule has when the system is in a given state. It follows that such maps are time varying, i.e., they in general specify different reactivities in correspondence of different temporal steps. Formally speaking, for each rule rwe define a reaction map Frthat maps states into real numbers: Fr:Q−→ R.(3) Since qis defined at any temporal step, the application of a reaction map Fr ultimately results in a positive real number that we will take as the reactivity of r in q. Such maps allow for a wide choice of possible definitions depending on the biological phenomenon under analysis: according to the traditional formulation of dynamic system it is not restrictive to consider real functions that in their structure include the state of the system plus factors such as the reactivity constants mentioned at the beginning of this section. As an example, consider a rewriting system having an alphabet made of five symbols, A={A, B, C, D, E}, and two rules, rand s: r:ABB kr →AC s:AE ks →BD (4) in which, consistently with notations traditionally adopted in biochemistry, we have specified constant reactivities (krand ks) that are peculiar to each rule— they could be, for example, kinetic parameters related to the chemical reactions respectively associated to these rules. Possible structures of the reaction maps might be the following ones: •simple reactivity constants Fr=kr Fs=ks •reactivities driven by the law of mass action Fr=krq(A)q(B) Fs=ksq(A)q(E) •reactivities depending on only the largest number of objects in the system that are visible to the rule Fr= max{q(A), q(B)} Fs= max{q(A), q(E)} •reactivities depending on an external promoter, like an enzyme capable of activating the reaction Fr=q(D) Fs={q(D)}2 Metabolic Algorithm with Time-varying Reaction Maps 47 In the following we will pick up this example as long as we need to illustrate the principles of the algorithm. 2.3 Reaction Weights Reaction maps are not used directly, as reactivities. Rather, their activity is proportionally distributed among the rules by means of so-called reaction weights. Every reaction weight then gives, for each symbol, a population amount a rule applies to in order to proportionally consume the corresponding object. By denoting with α(i) the ith symbol in a string α, with |α|the length of the same string, and with |α|Xthe number of occurrences of Xin α, then we define the reaction weight Wr¡αr(i)¢for r:αr→βrwith respect to the symbol αr(i). Normalization can be straightforwardly expressed in quantitative terms if we think that all rules co-operate, each one with its own reactivity, to consume all available objects. Thus, it must be: X ρ∈R|X∈αρ Wρ¡X¢= 1 ∀X∈ A (5) that is, for each symbol the sum of the reaction weights made over the rules containing that symbol in their left part equals unity. Holding this constraint, we can define the reaction weights for each r∈Ras Wr¡αr(i)¢=Fr X ρ∈R|αr(i)∈αρ Fρ , i = 1, . . . , |αr|(6) Note that, similarly to what happens in (5), we sum at the denominator over the rules containing the symbol αr(i) in their left part. Returning to our example, we have to compute Wr(A), Wr(B), Ws(A), Ws(E): Wr(A) = Fr Fr+FsWr(B) = Fr Fr= 1 Ws(A) = Fs Fr+FsWs(E) = Fs Fs= 1 (7) 2.4 Limitation, Rounding and State Transition For what we have said in the above, the available objects are consumed proportionally to the reaction weights. Then, in our example we have to choose whether to consume Wr(A)q(A) or Wr(B)q(B) (8) objects using r, and Ws(A)q(A) or Ws(E)q(E) (9) 48 L. Bianco, F. Fontana, V. Manca using s—provided for simplicity that all values in (8) and (9) are integer. The right choice is figured out by considering that every rule cannot consume more than the amount of the (reactant) object, taken with its own multiplicity in the reaction, whose availability in the system is lowest. Limitation, then, comes out for every rule by minimizing among all reactants participating to it: Λr= min i=1,...,|αr|nWr¡αr(i)¢q¡αr(i)¢ ¯¯αr¯¯αr(i)o(10) Still, Λris a real number. As opposite to this, a genuine object-based rewriting system must restrict the rule application domain to integer values. Instead of, for instance, rounding the minima obtained by (10), we prefer the following policy (later we will understand why): for every rule, compare the fractional part frac(Λr) of Λrto a random variable vrdefined between 0 and 1, and choose the floor of Λr if this fraction is smaller, the ceiling otherwise. In this way, new rounded minima result to be equal to: Λr=½floor(Λr),frac(Λr)≤vr ceil(Λr),frac(Λr)> vr .(11) As a result of this step we obtain the set {Λr, r ∈R}, containing the number of objects each rule will be applied to. In conclusion, for every symbol X∈ A the change in the number of objects due to ris equal to the stoichiometric factor of r, equal to |βr|X− |αr|X, times the value Λr: ∆r(X) = Λr(|βr|X− |αr|X) (12) It descends that for every symbol X∈ A the state evolves according to the following formula: qt+1(X) = qt(X) + X r∈R ∆r(X) (13) Again in our example, let us suppose that at time tit is q(A) = q(B) = q(C) = q(D) = q(E) = ˜q, furthermore Fr= 3/4Fs. Then, Λr= minn3/4Fs 3/4Fs+Fs ˜q, 1 2˜qo=3 7˜q Λs= minnFs 3/4Fs+Fs ˜q, ˜qo=4 7˜q After rounding Λrand Λs(here, for simplicity, we suppose to have found integers already at the limitation step) we have ∆r(A) = 3 7˜q(1 −1) = 0 ∆r(B) = 3 7˜q(0 −2) = −6 7˜q Metabolic Algorithm with Time-varying Reaction Maps 49 ∆r(C) = 3 7˜q(1 −0) = 3 7˜q ∆s(A) = 4 7˜q(0 −1) = −4 7˜q ∆s(B) = 4 7˜q(1 −0) = 4 7˜q ∆s(D) = 4 7˜q(1 −0) = 4 7˜q ∆s(E) = 4 7˜q(0 −1) = −4 7˜q in a way that qt+1(A) = ˜q−4 7˜q=3 7˜q qt+1(B) = ˜q+³−6 7+4 7´˜q=5 7˜q qt+1(C) = ˜q+3 7˜q=10 7˜q qt+1(D) = ˜q+4 7˜q=11 7˜q qt+1(E) = ˜q−4 7˜q=3 7˜q From the last equations it follows that X X∈A qt+1(X) = 32/7 ˜q < X X∈A qt(X) = 5˜q Interesting to see, in this system the total number of objects cannot increase along time. In other words in our example the following relation holds: X X∈A qt+1(X)≤X X∈A qt(X) 3 Flexibility of the Algorithm The proposed algorithm has two basic access points where parameters can be put into: the reaction maps, and the stochastic properties of vr. •Reaction maps can be defined with relative freedom, and even changed during the process according to the specific phenomenon under study. Their activity, in fact, is in any case normalized by the reaction weights. Occasionally some maps may result in null values: in this case reaction weights might arise in the form 0/0, and proper strategies must be put into action to handle them properly. Reaction maps, in conclusion, enable the fine control of the macroscopic, i.e., deterministic part of the process. 50 L. Bianco, F. Fontana, V. Manca •Conversely, the statistics of vrhas consequences on the system behavior that become as more important, as fewer objects are present in the system. In other words it influences the stochastic part of the process, i.e., its unpredictability in front of individual drifts from the average behavior. Although further research must be carried out to shift the metabolic algorithm closer to stochastic methods used in biochemistry [17], nevertheless the control of vralready allows to handle, at least to some extent, an interesting property of most discrete population dynamics, according to which the decision taken by an individual becomes as more crucial, as less populated the system is [6]. This feature is evident in the simulation of the Lotka-Volterra dynamics proposed in the following. The rounding policy expressed by (11) does not prevent that the resulting application of rules exceeds the available resources in the system. As an example suppose that, during a transition, it happens that Λr≥Λrfor each r∈R: in this case it is likely that the consequent application of the rules over-consumes at least some objects available in the system. To prevent this we must check that X r∈R Λr|αr|X≤q(X)∀X∈ A (14) otherwise the set of minima must be computed again. In first approximation vrcan be chosen to have a uniform distribution. We will present here an example in which a different choice of the random variable leads to more accurate simulation results. 3.1 Transparent Rules The metabolic algorithm allows to tune the activity of rules. Tuning is achieved by adding in the system so-called transparent rules, i.e., rules in the form αr→αr. Proper reaction maps can be selected to put such transparent rules in concurrence with the other, effective rules sharing common reactants. In this way, during a transition of the system every rule is applied as less intensively, as larger the reactivity value expressed by a concurrent transparent rule is. For instance, let us add a rule t:Akt →Ain the system expressed by (4). This leads to the following reaction weights: Wr(A) = Fr Fr+Fs+FtWr(B) = Fr Fr= 1 Ws(A) = Fs Fr+Fs+FtWs(E) = Fs Fs= 1 (15) Note that we can omit to compute Wt(A) due to the transparency of the corresponding rule. Clearly, Fttunes the action of rand sover the symbol A. In the limit case Ft=∞the rule tinhibits the action of rand sover A, since in this case we have Wr(A) = Ws(A) = 0. Finally, a rule t:ABE →ABE would inhibit both rand s in the same limit case. Metabolic Algorithm with Time-varying Reaction Maps 51 Transparent rules add further flexibility to the algorithm. In particular, they allow to observe the system evolution with the desired degree of resolution regardless of any consideration about the granularity of the temporal step. Changes in resolution are instead obtained by “hiding” objects to the system evolution by means of transparent rules. The way transparent rules work reflects an inherent attitude of the metabolic algorithm to scale its own resolution not along the time dimension, i.e., by means of a temporal scaling factor as it happens in most numerical methods. Rather, resolution is scaled by adapting the size of populations to the degree of precision expected for the experiment. 4 Results We show results coming from the predator-prey population dynamics, the Brusselator, and the PKC activation process. 4.1 Predator-Prey Population Dynamics The classic Lotka-Volterra population dynamics [10, 21] can be described by a simple set of rewriting rules in which Xare preys and Ypredators: r:Xkr −→ XX prey reproduction s:XY ks −→ Y Y predator reproduction t:Ykt −→ λpredator death (16) Here, we can tune the activity of every rule by selecting proper reactivity constants kr,ksand ktproportional to the rate of reproduction and death of both predators and preys. We postulate Fsto be constantly proportional to kstimes the maximum number between preys and predators, max{q(X), q(Y)}, that are present in the population at any system transition. Conversely, the remaining reaction maps are set to be constantly equal to the corresponding reactivity constants: Fr=kr Fs=ksmax{q(X), q(Y)}(17) Ft=kt Moreover we add transparent rules accounting for preys that are not reproducing or being consumed and for predators that are not eating or dyeing: u:Xku −→ Xprey standing by v:Ykv −→ Ypredator standing by (18) The set of reaction weights is then equal to 58 L. Bianco, F. Fontana, V. Manca 0 5 10 15 20 25 30 0 10 20 30 40 50 60 70 80 90 100 step concentration Ca PKC−a DAG PKC−i Fig. 4. PKC activation dynamics. The order of the elements in the legend is the same as the order of their final concentrations within the plot. meanwhile we associate a constant reactivity map to each transparent rule, in our case F= 50. These reactivity maps are quite simple but in the future we intend to investigate the effectiveness of more complex reactivity maps in the case of the PKC model. The description of the whole set of weights Wri, that can be calculated object by object in the way introduced in previous sections, is omitted. Rather, we present some simulation results obtained using our algorithm. In Figure 4 we see that, in accordance with results obtained in [2], PKC −idecreases to zero while PKC −agrows up until reaching a stationary maximum. Figure 5 represents the characteristic dynamics of the diacylglycerol-protein kinase C (DAG.PKC) complex. 5 Discussion All rules discussed so far do not present any target specification. This aspect needs further discussion due to its importance. Let’s consider the following rule r, present in a membrane wi: Metabolic Algorithm with Time-varying Reaction Maps 59 0 10 20 30 40 50 60 70 80 90 100 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 step concentration DAGPKC Fig. 5. DAG.PKC complex dynamics. AB →BINjC , Fr where Fris the reactivity map associated to r. Its meaning is the following: whenever Ajoins Binside wi, they combine and produce an object Cinside the same membrane, meanwhile an object Bleaves wiand reaches the membrane wj. In such a way raffects objects that are present in two different membranes. In particular, from a structural viewpoint, the elements Bthat are present in wihave to be distinguished from the elements Bthat are present in wj(and, in fact, this is the effect of compartmentalization). For this reason roriginates four metabolic equations describing the behavior of its four distinct elements: ∆r(Awi) = −Λr,wi ∆r(Bwi) = −Λr,wi ∆r(Bwj)=+Λr,wi ∆r(Cwi) = +Λr,wi where we have introduced the label of the membrane containing every element as subscript. In this way we can see that the variation due to ron the objects Bplaced inside wj, i.e., ∆r(Bwj), depends on the concentrations of Aand Blocated in wi 60 L. Bianco, F. Fontana, V. Manca as stated by the subscript notation (note that this dependence is hidden behind the Λrfactor). This simple evolution rule is powerful enough to show that movements of objects between membranes can be handled easily by considering, as distinct elements, two objects of the same type located in different regions. This additional information introduces a notational overhead of targeting every object with the label of the membrane containing the respective object. On the other hand it does not introduce any conceptual complication. The case in which elements appearing in the antecedent of the rule are placed inside different membranes can be handled similarly. The only difference that must be taken into account is that the set of weights has to be calculated by considering, in principle, the whole set of rules of the P system rather than the set of rules of a single membrane. 6 Conclusion and Ongoing Research Systems biology demands for novel procedures capable of representing biological processes with both accuracy and flexibility. In front of a huge and well-rooted family of numerical schemes, traditionally devoted to figure out the dynamics of systems described by differential equations, alternative algorithms based on a symbolic representation of the phenomena promise to deal more naturally with the structural characteristics of the biomolecules and with the biochemical reactions such molecules give rise to. By using the same kind of representation, our algorithm moreover seems to handle in a straightforward and efficient way those conditions in which few molecules have an important impact on the system dynamics, where most traditional numerical strategies are no longer effective and must be substituted by stochastic algorithms. Successful simulations conducted on two paradigmatic nonlinear processes in biochemistry, namely the Lotka-Volterra population dynamics and the BZ reaction, plus the experiment conducted with the PKC activation process, ask for further test the potential of the metabolic algorithm. Our present research aims to simulate some fundamental signal transduction networks, in particular the PER and TIM cycle in the circadian oscillation in Drosophila. References 1. D. Besozzi, G. Ciobanu: A P system description of the sodium-potassium pump. In Membrane Computing, 5th International Workshop, WMC 2004 (G. Mauri, G. P˘aun, M.J. P´erez–Jim´enez, G. Rozenberg, A. Salomaa, eds.), LNCS 3365, Springer-Verlag, Berlin, 2005, 210–223. 2. U.S. Bhalla, R. Iyengar: Emergent properties of networks of biological signaling pathways. Science, 283 (January 1999), 381–387. 3. L. Bianco, F. Fontana, G. Franco, V. Manca: P systems for biological dynamics. In [4]. Metabolic Algorithm with Time-varying Reaction Maps 61 4. G. Ciobanu, Gh. P˘aun, M.J. P´erez–Jim´enez, eds.: Applications of Membrane Computing. Springer, Berlin, 2005. 5. I.R. Epstein, K. Showalter: Nonlinear chemical dynamics: Oscillations, patterns, and chaos. J. Phys. Chem., 100, 31 (1996), 13132–13147. 6. A. Goldbeter: Computational approaches to cellular rhythms. Nature, 420 (November 2002), 238–244. 7. H. Kitano: Computational systems biology. Nature, 420 (November 2002), 206–210. 8. J.C. Leloup, A. Goldbeter: A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. J. of Biological Rhythms, 13, 1 (1998), 70–87. 9. A. Lindenmayer: Mathematical models for cellular interaction in development. Part I and part II. J. of Theoretical Biology, 18 (1968), 280–315. 10. A.J. Lotka. Undamped oscillations derived from the law of mass action. J. Am. Chem. Soc., 42 (1920), 1595–1599. 11. V. Manca, L. Bianco, F. Fontana: Evolutions and oscillations of P systems: Theoretical considerations and applications to biochemical phenomena. In Membrane Computing, 5th International Workshop, WMC 2004 (G. Mauri, Gh. P˘aun, M.J. P´erezJim´enez, G. Rozenberg, A. Salomaa, eds.), LNCS 3365, Springer-Verlag, Berlin, 2005, 63–84. 12. A.C. Newton: Protein Kinase C: Structure, function, and regulation. J. Biol. Chem., 270 (1995), 28495–29498. 13. K. Ohkusu: Elucidation of the protein kinase C-dependent apoptosis pathway in distinct subsets of T lymphocytes in mrl-lpr/lpr mice. Eur. J. Immunol., 25 (995), 3180. 14. M.J. Perez-Jimenez, F.J. Romero-Campero: Modelling EGFR signalling network using continuous membrane systems. Third Workshop on Computational Methods in Systems Biology, Edinburgh, 2005. 15. P. Prusinkiewicz, M. Hammel, J. Hanan, R. Mech: Visual models of plant development. In Handbook of Formal Languages (G. Rozenberg, A. Salomaa, eds.), volume III: Beyond Words, Springer-Verlag, Berlin, 1997, 535–597. 16. Gh. P˘aun: Computing with membranes. J. Comput. System Sci., 61, 1 (2000), 108– 143. 17. C.V. Rao, D.N. Wolf, A.P. Arkin: Control, exploitation and tolerance of intracellular noise. Nature, 420 (November 2002), 231–237. 18. G. Rozenberg, A. Salomaa, eds.: Handbook of Formal Languages. Springer-Verlag, Berlin, 1997. 19. I. Stamatopoulou, M. Gheorghe, P. Kefalas: Modelling dynamic organization of biology-inspired multi-agent systems with communicating x-machines and population p systems. In Membrane Computing, 5th International Workshop, WMC 2004 (G. Mauri, Gh. P˘aun, M.J. P´erez-Jim´enez, G. Rozenberg, A. Salomaa, eds.), LNCS 3365, Springer-Verlag, Berlin, 2005, 389–403. 20. Y. Suzuki, H. Tanaka: Chemical oscillation in symbolic chemical systems and its behavioral pattern. In Proc. International Conference on Complex Systems (Y. BarYam, ed.), Nashua, NH, September 1997. 21. V. Volterra: Fluctuations in the abundance of a species considered mathematically. Nature, 118 (1926), 558–560.