# Metabolic P systems

Post-publication activity

Curator: Vincenzo Manca

Metabolism is one of the basic phenomenon on which life is based. Any living organism has to maintain processes which: i) introduce matter of some kinds from the external environment, ii) transform internal matter by changing the molecule distribution of a number of biochemical species (substances, metabolites), and iii) expel outside matter which is not useful or dangerous for the organism. Molecule distribution identifies the metabolic state of the system in question, and can be represented as multiset of type $$n_A A + n_B B + \ldots n_Z Z$$ giving the numbers $$n_A, n_B, \ldots , n_Z$$ of molecules for each of molecular species $$A, B, \ldots , Z\ .$$ Of course, life cannot be reduced to metabolism, but no life can exist without such a kind of basic mechanism. Metabolic P systems or P metabolic systems (shortly MP systems) represent metabolic processes in a discrete mathematical framework. A metabolic P system is essentially a multiset grammar where multiset transformations are regulated by functions. Namely, a multiset rule like $$A+B \to C$$ means that a number $$u$$ of molecules of kind $$A$$ and the same number $$u$$ of molecules $$B$$ are replaced by $$u$$ molecules of type $$C\ .$$ The value of $$u$$ is the flux of the rule application. Assume to consider a system at some time steps $$0, 1,2, \ldots, i\ ,$$ and consider a substance $$x$$ that is produced by rules $$r_1, r_3$$ and is consumed by rule $$r_2\ .$$ If $$u_1[i], u_2[i], u_3[i]$$ are the fluxes of the rules $$r_1, r_2, r_3$$ respectively, in the passage from step $$i$$ to step $$i+1\ ,$$ then the variation of substance $$x$$ is given by$x[i+1] - x[i] = u_1[i] - u_2[i] + u_3[i]\ .$

Applying this method of mass partition among rules, for all the substances, we can compute the dynamics of a system from the fluxes of rules at every step. In an MP system, it is assumed that, in any state, the flux of each rule is provided by a function, called regulator of the rule. The letter P of MP systems comes from the theoretical framework of P systems introduced by Gheorghe Paun (2002), in the context of membrane computing. In fact, MP systems are a special class of P systems introduced in 2004 (see V. Manca, L. Bianco, F. Fontana (2005)) for expressing metabolism in a discrete mathematical setting.

## Multiset rewriting: vector notation and molar perspective

P systems, which provide the right conceptual framework for MP systems development, were introduced in the context of Formal Language Theory for defining a computation model inspired by biological cells (in the same research line of Lindenmayer's L systems and Head's H systems). A crucial aspect of P systems is the passage from the classical notion of string rewriting rule, to that one of multiset rewriting rule. This new perspective changes completely the typical mathematical setting of languages, grammars and automata, based on strings, and opens a new perspective where rewriting rules become a natural representation of (bio)chemical reactions ('rule' and 'reaction' will be used synonymously in the following). If a rule transforms a multiset of symbols/molecules, then when a conventional order, for example the alphabetic order, is fixed over a set $$S$$ of $$n$$ substances (molecule types), then a rule such as $$2A+B \to C$$ is easily represented by two vectors of $$N^n\ ,$$ that is, a left vector $$(2, 1, 0, \ldots , 0)\ ,$$ and a right vector $$(0, 0, 1, \ldots , 0)\ .$$ In general for a rule $$r$$ its left vector is denoted by $$r^-$$ and its right vector is denoted by $$r^+ \ .$$ The vector difference $$r^\# = r^+ - r^-$$ is the stoichiometry balance of rule $$r\ .$$ This vector reading of rules is only a notational change and does not alter, in any way, the notion of rule of P systems (that we will continue to use), nevertheless it is very useful for applying basic concept of linear algebra to the analysis of MP systems. Another important aspect of multiset rewriting is its natural reading according to the usual chemistry use. In fact, if we fix a number of population size, called conventional mole, for example of $$10000$$ units, then a rule $$2A+B \to C\ ,$$ when it is applied with a reaction flux of $$3.2$$ moles, it transforms $$64000$$ molecules $$A$$ and $$32000$$ molecules $$B$$ into $$32000$$ molecules $$C\ .$$ Figure 1: An MP graph. Nodes: triangles represent matter introduction and expulsion, big circles stand for substances, small full circles for reactions labeled by mathematical expressions denoting regulators. Edges: transformation edges go from substances to reactions (consumption) and from reactions to substances (production), regulation red edges go from substances to reactions and represent the substances on which each regulator depends (also called tuners of that regulator).

## MP systems: a formal definition

Substances, reactions, and regulators are essential components of an MP system. Parameters are added to these components, when magnitudes different from substances are essential in regulating reactions (the term magnitude will refer to a substance or to a parameter). Moreover, a temporal interval $$\tau$$ and a conventional mole size $$\nu$$ are considered, which specify the time and population (discrete) granularities respectively. Finally, for each substance, its mole mass is specified.

An MP system is a discrete dynamical system where dynamics is considered at steps indexed in the set $$N$$ of natural numbers. The temporal interval, the mole size, and substance mole masses belong to the set $$R$$ of real numbers. They are scale factors and do not enter directly in defining the dynamics of a system, but are essential for interpreting it at a specific physical level of mass and time granularity.

From a mathematical point of view, an MP system $$M$$ of type $$(n,m,k)\ ,$$ that is, of $$n$$ substances, $$m$$ reactions and $$k$$ parameters (this type will be implicitly assumed), is specified by (see also V. Manca (2010)):

• a set $$S$$ of substances, considered in a conventional order, determining, for any metabolic state of the system, a vector $$X$$ of $$R^n$$ with components expressing substance quantities in conventional mole units;
• a set $$R$$ of reactions, which are given by $$m$$ pairs $$(r^-_1,r^+_1) \ldots , (r^-_m, r^+_m)$$ $$\in N^n \times N^n \ ,$$ composed by the left and right vectors of the reactions (relative to the reactants and to the products respectively). The matrix $$A = (r^\#_1, \ldots , r^\#_m)$$ is the stoichiometric matrix associated to the reactions having as columns the stoichiometry balances of the rules;
• a function $$H: N \to R^k$$ providing, at each step $$i \in N\ ,$$ the vector $$H[i]$$ of parameters;
• a vector of regulators (or flux regulation functions) $$\Phi = (\varphi_1, \ldots , \varphi_m)\ ,$$ where $$\Phi : R^{n} \times R^{k} \to R^m$$ provides the fluxes of reactions corresponding to any global state of the system, that is, a pair in $$R^n \times R^k$$ constituted by the metabolic state and by the parameter vector. However, given a reaction $$r\ ,$$ only some of the substances or parameters, which are called the tuners of the reaction, may occur as arguments of the corresponding regulator $$\varphi_r\ .$$ Some constraints may be imposed to the fluxes provided by regulators, which may be of general nature, or may be specific to some classes of systems (for example, fluxes should not be negative, and the sum of fluxes of all reactions consuming a substance $$x$$ cannot exceed the quantity of $$x$$).
• a time interval $$\tau \in R$$ between two consecutive steps;
• a conventional mole $$\nu \in R\ ;$$
• a vector $$\mu \in R^n$$ of the mole masses of substances.

Given a vector $$X \in R^n$$ relative to an initial state of a given system, the dynamics of $$M$$ is specified by the following vector recurrent equation, called $$EMA[i]$$ (Equational Metabolic Algorithm), where $$\times$$ is the usual matrix product, and, in dependence on the context, $$+$$ is the usual sum or the component-wise vector sum$\tag{1} X[i+1] = A \times U[i] + X[i]$

providing the state of the system $$X[i+1]\ ,$$ for each step $$i \in N\ ,$$ by means of the vector of fluxes $$U[i]= (u_r[i] \;|\; r \in R)$$ where $$(u_r[i] = \varphi_r(a[i], b[i], \ldots )$$ and $$a[i], b[i], \ldots$$ are components of $$X[i], H[i]$$ which are the tuners of reaction $$r\ .$$

MP dynamics, computed by EMA, is nothing else than the difference equation method, suitably devised in the context of discrete metabolic transformations.

## MP graphs and MP grammars

The intuitive meaning of an MP system can be represented by an MP graph (see figure 1, V. Manca, L. Bianco (2008)) expressing a hydrodynamic analogy: substance nodes are vessels of different liquids and transformation edges are pipes connected by means of reaction nodes (liquid passages imply also liquid transformations) where regulators act as taps which regulate the flux of liquids passing through the pipes, in dependence on the values assumed by their tuners. Apart the scale factors, all the elements specifying an MP system are represented in its MP graph (parameter evolution functions can be add as labels of parameter nodes). Moreover, the whole information inside an MP graph can be also expressed by an MP grammar, that is, a set of multiset rewriting rules, each of them equipped by a corresponding regulator (and parameters equipped by their evolution functions). The following is an example of MP grammar, where $$\emptyset$$ denotes an empty multiset ($$\emptyset$$-rules correspond to introduction or expulsion of matter) and substance symbols occurring in regulators denote the corresponding substance quantities.

• $$r_1 : \emptyset \to A \;\;\;\; \varphi_1 = 0.5 p$$
• $$r_2 : \emptyset \to B \;\;\;\; \varphi_2 = 0.2 C$$
• $$r_3 : B \to \emptyset \;\;\;\; \varphi_3 = 0.1$$
• $$r_4 : A \to C \;\;\;\; \varphi_4 = 0.6 A/C$$
• $$r_5 : C \to\emptyset \;\;\;\; \varphi_5 = 0.4$$
• $$p= 0.2 \;\;\;\; p[i+1]= p[i]+ 0.2$$

MP graph notation, is similar to typical graphical notations adopted in many biochemical and biological contexts. Nevertheless, in MP graphs, function labels for regulators can express, in a uniform and rigorous way, very general relations of dependence among reactions, substances and parameters. Beside their direct reading as graphs and grammars, MP systems have be shown to be equivalent to a class of continuous Petri nets (see A. Castellini, G. Franco, V. Manca (2009)). However, the main specific interest of MP systems, with respect to other related formalisms, is their natural possibility of developing theories and tools for deducing MP models from time series of given metabolic systems, as it will be shown in the next section.

In the given definition of MP system, the classical compartmentalization of P systems is missing, but it could be easily added, in many ways, with a further complication of the structure. A simple way for dealing with a (static) membrane structure is that of considering the same substance $$x\ ,$$ possibly occurring in two different membranes $$L, M\ ,$$ as two distinct substances $$x_L\ ,$$ and $$x_M\ .$$ In this way, the passage of $$x$$ from $$L$$ to $$M$$ is modeled as a transformation of $$x_L$$ into $$x_M$$ (see V. Manca (2010) for more details on this membrane representation). In this perspective, the local character of a reaction inside a membrane is naturally translated by the dependence of its regulator only from substances and parameters which are internal to the membrane. According to this representation, the theory of log-gain, which will be outlined in the next section, can be entirely applied even when a membrane structure is relevant in the metabolic system under investigation.

Here, rather than on substance localizations, we focus on the specific aspect of quantitative dynamics. However, at a deeper analysis, a particular notion of membrane turns to be crucial in the MP systems too. In fact, in an MP graph nodes can be seen as compartments where matter flows by changing its biochemical localization instead of its spatial localization. In this sense MP systems resemble another model, defined in the context of membrane computing, where spikes move between neurons, interconnected as nodes of a graph, and the role of objects is played by the signals of communicating neurons (especially in the non-synchronized neural spiking systems, see M. Cavaliere, O. Ibarra, G. Paun, O. Egecioglu, M. Ionescu, S. Woodworth (2009)).

## Inverse dynamical problems and log-gain theory

Given an MP system, the recurrent equation system $$EMA$$ generates the evolution of substances, according to the MP grammar of the system (starting from an initial metabolic state, and with the knowledge of parameter evolution). In other words, as it is illustrated by figure 1, when regulators $$\Phi$$ are given, then fluxes can be computed and then the substance variations follow easily from the stoichiometry. The inverse dynamical problem is the opposite verse of this process. In fact, let us assume to 'observe' a metabolic system for a number of steps (separated by some temporal interval). How can we discover an MP system which provides the dynamics which we observe? The substance and the reactions in question are given by the particular phenomenon we want to describe. Moreover, the stoichiometry can be deduced by a basic knowledge of the phenomenon, but how to know the fluxes of matter transformation which are responsible of the observed evolution? This is the problem of flux discovery.

### Flux discovery

If we know the temporal series of substance quantities, that is the vector sequence $$(X[i] | i\in N)\ ,$$ then we can read the systems $$(EMA[i] | i \in N)$$ by reversing the known values with the unknown ones. In fact, by denoting the difference $$X[i+1] - X[i]$$ by $$\Delta X[i]$$ and assuming $$n$$ substances and $$m$$ reactions, then we get the following system $$ADA[i]$$ (Avogadro and Dalton Action) of $$n$$ equations and $$m$$ unknowns (the $$m$$ components of the flux vector $$U[i]$$)$\tag{2} A \times U[i] = \Delta X[i]\ .$

The acronym $$ADA$$ is used because, following the Avogadro principle of chemistry, in any reaction the same flux of matter involves all the substance of the reaction according to the stoichiometric coefficients of the reaction; moreover the variation of any substance is obtained by means of the algebraic sum all the fluxes of reactions consuming and producing that substance, that is according a law recalling Dalton principle in chemistry (see also V. Manca (2008a)).

However, the equations (2) above are not enough for solving the flux discovery problem. In fact, usually the number of reactions is greater than the number of substances, therefore the $$ADA$$ systems are not univocally solvable. For obtaining such a solvability we need to add more equations to $$ADA\ .$$ The log-gain principle satisfies this necessity by postulating, in terms of equations, that the relative variation of a (strictly positive) reaction flux has to be a linear combination of the relative variations of the tuners of the reaction. This statement corresponds to a special formulation of the general biological principle of allometry (see L. Bertanlaffy (1967)), which claims a harmonic form in biological changes. In technical terms, given a time series $$(Z[i] \;|\; i \in N)\ ,$$ its (discrete) log-gain $$\Gamma(Z(i))$$ is given by$\Gamma(Z(i)) = \frac{Z(i+1) - Z(i)}{Z(i)}$

(if $$Z(i)$$ is a vector, difference and division operations are intended to be applied in component-wise manner).

Therefore, the log-gain equations can be written in the following manner, where now $$(Z[i]\;|\; i \in N)$$ is the time series of magnitudes of a given MP system (of type $$(n, m, k)$$), $$U[i] = (u_r[i] \; |\; r \in R)\ ,$$ and for each $$r \in R\ ,$$ row vector $$B_r[i]$$ is given by the tuner log-gain coefficients of $$r$$ (coefficients of $$B_r[i]$$ are null for magnitudes which are not tuners of $$r$$)$\tag{3} \Gamma(u_r[i]) = B_r\times \Gamma(Z[i])$

of course, if $$B$$ is the matrix $$m \times n+k$$ of vectors $$B_r \ ,$$ then$\tag{4} \Gamma(U[i]) = B \times \Gamma(Z[i])\ .$

The term 'log-gain' is used because, in differential context, the infinitesimal relative variation of a magnitude corresponds to the derivative of its logarithm. Elaborating on this principle it can be shown (see V. Manca (2007), V. Manca (2008), V. Manca (2009), V. Manca (2009a)) that a system of equations, can be obtained by adding log-gain equations to $$ADA[i]\ .$$ Under very wide hypotheses, the resulting system turns out to be univocally solvable (the main hypothesis involves the rank of the stoichiometric matrix). An essential aspect of the resulting system is a kind of step displacement of the log-gain equations with respect to the $$ADA$$ equations. In fact, the latter are the equations of $$ADA[i+1]\ ,$$ while the former relate to equations (3). This displacement implies an inductive solution but avoids nonlinearity. We only mention here three essential sub-problems related to the flux discovery problem, which were solved in all the specific cases we considered, but for which general resolution strategies are under investigation: i) the determination of the initial flux vector $$U\ ,$$ ii) the determination of the tuners of each reaction, and iii) the determination of suitable linearly independent vectors among the stoichiometric balances of the $$m$$ reactions (see G. Franco, V. Manca, R. Pagliarini (2009), R. Pagliarini, G. Franco, V. Manca (2009)).

### Regulator discovery

Now, let us assume that we discovered the reaction fluxes at each step, that is, for each $$i >0$$ we got the flux vector $$U[i]\ .$$ We have to solve a final problem in order to complete the determination of an MP system underlying the observed dynamics, that is, an MP system able to generate the substance time series we observed. This last problem is the regulator discovery, that is, the determination of some functions $$\Phi$$ providing the fluxes we got from the solution of the flux discovering problem. At this end, we used regression techniques such as Least Square Method, Stepwise Regression, Artificial Neural Networks (and integrations and specializations of them, see G. Franco, V. Manca, R. Pagliarini (2009), A. Castellini, V. Manca (2009)). For example, assume to fix the forms of functions which could approximate regulators, say as polynomials of a certain type $$P_r(A, B, \ldots, c_1, c_2, \ldots, c_{kr})$$ (of tuners $$A, B, \ldots$$ and coefficients $$c_1, c_2, \ldots, c_{kr}$$). This means that we need to find $$k1$$ coefficients for the first regulator, $$k2$$ coefficients for the second regulator, and so on. Therefore, by using a number $$T$$ of observation steps greater than $$k1+k2+ \ldots km\ ,$$ the approximation of the right regulators can be reduced to the least square approximation problem associated to the following system of equation, for $$1 \leq i \leq T\$$ and $$r \in R\:$$

$$\tag{5} u_r[i] = P_r(A[i], B[i], \ldots, c_1, c_2, \ldots, c_{kr}) .$$

### Log-gain validation

The methodology outlined in the previous section was validated in the following way. We started from an MP system $$M$$ with regulators $$\Phi\ ,$$ defined according to the definition of section MP systems: a formal definition, then we generated the time series $$(X[i] \;|\; i>0 \in N)$$ of $$M$$ from a chosen initial state $$X$$ by using the $$EMA$$ system. At this point, we applied the log-gain method for discovering the flux time series $$(U[i] \;|\; i>0 \in N)\ ,$$ and then we solved the regulator discovering problem for these fluxes, getting regulators $$\Psi\ .$$ The discovered regulators $$\Psi$$ resulted to be a good approximation of $$\Phi\ .$$ The whole validation process was successful in all the cases we tested, even if with different, but always reasonable, approximations (see L. Bianco, F. Fontana, G. Franco, V. Manca (2006), V. Manca, R. Pagliarini, S. Zorzan (2009), A. Castellini, V. Manca (2009), Official site of MetaPlab software), and Log-gain principles resulted to be essential for devising specific regression techniques based on Stepwise Regression, and Artificial Neural Networks. Log-gain method was also applied in the approximation of periodical real functions, by reconstructing, from a number of time series, MP grammars which provided very complex behaviours with a surprising high level of precision (see V. Manca, L. Marchetti (2010)).

## MP biological modeling

P systems are a natural tool for biological modeling, especially of qualitative nature (see for example G. Ciobanu, G. Paun, M.J. Perez-Jimenez eds. (2006) and G. Franco, V. Manca (2003)). MP systems emphasise a perspective of discrete quantitative modeling. In fact, quantitative dynamical models of important phenomena were reconstructed in terms of MP systems, e. g., Goldbeter’s mitotic oscillator (see figure 2), Belousov-Zhabotinski reaction in the Nicolis and Prigogine’s formulation, Lotka-Volterra’s Prey-Predator model (see L. Bianco, F. Fontana, G. Franco, V. Manca (2006)), and NPQ (Non-Photochemical Quenching in photosynthesis, see V. Manca, R. Pagliarini, S. Zorzan (2009)). In all these cases a complete concordance with the classical models was found. Moreover, some synthetic oscillators with interesting behaviors were easily discovered (see section MP dynamics and MP oscillators, and MP models for NonPhotochemical Quenching (where no standard mathematical model is known) was directly deduced, by experimental data, by using the Log-gain theory (see V. Manca, R. Pagliarini, S. Zorzan (2009)). Many MP models are available in the Official site of MetaPlab software (e.g., Mitoticus, Sirius, Brusselatus, Volteranus, Photosyntheticus).

## MP systems and ordinary differential equations

In a first approximation MP systems can be viewed as a discrete approximations of metabolic models based on ordinary differential equations ($$ODE$$), by means of recurrent difference equations. In fact, if we consider a very small temporal interval $$\tau$$ of an MP system, then fluxes approximate to differentials and $$EMA$$ reduces to a numerical integration performed according to the Euler method (see F. Fontana, V. Manca (2007)). However, when the time interval cannot be approximated to an infinitesimal time, MP models introduce a perspective which is radically different from the ODE perspective: it is essentially observational and more related to the transformations of the systems along the observation steps, rather than to the internal kinetic of reactions. Therefore, the MP approach is systemic vs kinetic, discrete vs continuous, global vs local. This different perspective seems to have computational and modeling advantages, and the notion of flux on which it is based plays a completely different role with respect to the other classical approaches of systems biology (see E. Klipp, C. Wierling, A. Kowald, H. Lehrach, H. Herwig (2009)), by opening the possibility of algorithmic and algebraic procedures for discovering the regulation logic of metabolic systems.

In fact, in ODE models, and in stochastic models (see D.T. Gillespie (1977), D.T. Gillespie (2007)), a critical point is constituted by the evaluation of kinetic rates or stochastic parameters. This is a hard and unreliable task, because very often the exact microscopic dynamics of reactions is not completely known, or even when it can be completely studied in vitro, its in vivo behaviour could significantly differ from its experimental evidence. Moreover, differential laws are based on instantaneous effects (like in gravitational or electromagnetic fields), which, in biochemical contexts, is an implicit assumption that is satisfied when chemicals are supposed to be homogeneously distributed in the volume, or immediately at disposal to react. But, this assumption has no evident correspondence in many cases of biological interactions. On the other side, in MP models, the dynamical perspective is focused on the transformations between observation steps, rather than on reaction kinetics, and passages, from a state to another one, are considered at macroscopic time scales. Regulators distribute reagents among transformations (reactions), and if they work correctly (with respect to the reality they model) the quantities they assign are computed without no direct knowledge or assumption about the physical mechanism underlying the single biochemical processes. In this sense, no specific hypothesis of chemo-physical nature (about homogeneity, or concentration) is required for solving the inverse dynamics problem and for providing, via the log-gain method, the mathematical form of regulators. If they are adequate, then they express matter transformation rules (associated to a stoichiometry) which synthetise the systemic logic of the observed behaviour they model.

In other words, the MP approach assumes that some species are related by means of some transformations, and that the variations are due to the global action of these transformations. Their execution could be realized in many ways and could involve underlying very complex sub-transformations. However, this is outside the objective of the model. It only tries to explain what is observed in terms of the chosen species and transformations. If the choices are not the right one, this means that the model was not adequately defined, but this is independent from the methodology, it is only a matter of modeling design. In conclusion, MP modeling, is deliberately at a different, more abstract, level with respect to ODE models. This does not means that it is less adherent to the reality, but simply that it is focused on a different level of reality. The Flux discovery shows as the observational macroscopic approach to the inverse dynamical problem circumvents the critical point of kinetic rate evaluation of ODE. Of course, in this way many important dynamical details can be lost, but in many cases what is relevant is the global pattern of a behaviour and the main correlations among the involved magnitudes, and often a higher level of analysis is the only possible for deducing the logic underlying a phenomenon. Some variants of P systems, such as stochastic and probabilistic P systems, used in the modeling of cellular systems, are discrete approaches which may complement the MP perspective with a nondeterministic dynamical perspective (see M. Gheorghe, V. Manca, F.J. Romero-Campero (2009), F.J. Romero-Campero, M. Gheorghe, G. Ciobanu, J.M. Auld, M.J. Perez-Jimenez (2007), D. Pescini, D. Besozzi, G. Mauri, C. Zandron (2006)).

However, an important synergy can be established between ODE and MP models. In fact, let us suppose, to have an ODE model of a metabolic process. According to it, any derivative of substance quantity is the sum of some additive terms coming from the infinitesimal fluxes of the reactions. Now, if we consider a time interval $$\tau$$ and perform a numerical integration at a discrete time $$\Delta t\ ,$$ then the integration of the infinitesimal fluxes along $$\tau/\Delta t$$ steps ($$\tau$$ multiple of $$\Delta t$$), provides the MP fluxes of the reactions involved in the system in the time interval $$\tau\ .$$ This means that we get exactly what the log-gain theory provides by its flux discovery procedure. In other words, we get the macroscopic fluxes from the ODE microscopic ones. From these fluxes, by approximation techniques (Least Square Method, Stepwise regression, ANN) we can derive the flux regulation functions of an MP system which provides the same dynamics (along the steps at time interval $$\tau$$). It would be possible, that at this different temporal grain, some systemic effects emerge which could shed new light on the analysis of the modeled phenomenon.

## MetaPlab software

A Java software, called MetaPlab, was developed starting from a prototypal versions (see L. Bianco, V. Manca, L. Marchetti, M. Petterlini (2007)). MetaPlab is downloadable from the Official site of MetaPlab software. This platform enables the user to design MP models by means of some useful graphical tools, to simulate their dynamics, and to automatise the MP flux discovery and regulator discovery procedures. MetaPlab is based on an extensible set of plugins, namely Java tools, for solving specific tasks relevant in the framework of MP systems. Besides the basic plugins (dynamics computation, log-gain flux discovery, and linear regression) there are other plugins for exporting MP models in different formats (HTML, SBML, see V. Manca, L. Marchetti (2009)), for dynamics visualizations, for integrating differential equations, for computing from ODE integrations the MP fluxes of MP models, and for finding reaction tuners by means of ANN (see A. Castellini, V. Manca (2009)). A guide for this software V. Manca et. al. (2009) is available at Official site of MetaPlab software. Other plugins are under development: importation of biological networks from on-line databases, advanced regression tools, and determination of initial fluxes by specific optimization methods (see R. Pagliarini, G. Franco, V. Manca (2009)).

## MP dynamics and MP oscillators

MP systems are dynamical systems which seem flexible enough for revisiting in special context of biological relevance many general concepts of dynamical systems. In particular, notions such as of rigid and soft attractor, creods, and bistability can receive a natural description in terms of fluctuations and quasi-dynamics, where parameter driven regulations play an essential role (see V. Manca (2009b). MP systems are a natural framework for studying biological oscillators, which represent a very important class of metabolic dynamics. Many examples of MP oscillators were found (Sirius, Sirius geminus, Sirius ternarius, Mizar, Vega, see V. Manca (2009b), V. Manca (2009c)). Moreover, several metabolic schemata of biological oscillatory patterns can be easily represented by MP graphs. Recently, the MP system Goniometricus was found, having 3 rules and only linear regulators, where two parameters provide, by means of their dynamics, the sine and cosine functions with a precision of 14 decimal digits (see V. Manca, L. Marchetti (2010)). Figure 3: The dynamics of Sirius quaternarius with initial state$A= 100, B= 70, C= 100, D= 20\ .$

The following MP grammar identifies an MP system, called Sirius quaternarius, providing a perfect oscillatory pattern depicted in figure 3.

• $$r_1 : A + C \to 2A + D$$ $$\;\;\;\; \varphi_1 = 20A/p$$
• $$r_2 : A \to B$$ $$\;\;\;\; \;\;\;\; \;\;\;\;\;\;\;\; \varphi_2 = 0.14 B/p$$
• $$r_3 : B \to \emptyset$$ $$\;\;\;\; \;\;\;\; \;\;\;\;\;\;\;\; \varphi_3 = B/71$$
• $$r_4 : C \to \emptyset$$ $$\;\;\;\; \;\;\;\; \;\;\;\;\;\;\;\; \varphi_3 = 1.5 C/21.5$$
• $$r_5 : \emptyset \to C$$ $$\;\;\;\; \;\;\;\; \;\;\;\;\;\;\;\; \varphi_4 = 20$$
• $$r_6 : D \to \emptyset$$ $$\;\;\;\; \;\;\;\; \;\;\;\;\;\;\;\; \varphi_5 = 0.05$$
• $$p= 120 + 0.14 B\ .$$

## Extensions and applications of MP systems

Further research on metabolic P systems ranges from mathematical aspects to applications to themes of systems biology, computational synthetic biology, and complex systems analysis. On the mathematical side, we mention the possibility using Galois connections for a better understanding of substance-reaction interactions. In fact, given a substance $$x\ ,$$ we denote by $$R(x)$$ the set of reactions where $$x$$ occurs (as product or reactant), but symmetrically, given a reaction $$r\ ,$$ we can define $$S(r)$$ as the set of substances involved in the reaction r. If we extend $$R$$ and $$S$$ as functions from set of substances to set of reactions, and vice versa, we get a Galois connection, which is a very general and powerful algebraic concept. It seems possible that many metabolic concepts, are related to properties which can be analyzed in this algebraic setting. On the applicative side, we argue that many problems on protocells constitution and evolution, and the phenomenon of their autopoiesis are related to a sort of second-order dynamics which pushes metabolic systems to evolve. This process, and its relationship with the replication process of proto-biomolecular complexes is one of the possible main streams of research lines about MP systems. Finally, any process where some resources transform into other ones, and an income-outcome activity is kept for some internal finalities, falls into a typology which can be read sub specie metabolica. In this sense, applications to economic and financial activities could enlarge the modeling possibility of MP systems.