Inverse Dynamical Problems: An Algebraic Formulation Via MP Grammars

Inverse Dynamical Problems: An Algebraic Formulation Via MP Grammars Vincenzo Manca and Luca Marchetti University of Verona, Department of Computer Science Strada Le Grazie 15, 37134 Verona, Italy [email protected] [email protected] Summary. Metabolic P grammars are a particular class of multiset rewriting grammars introduced in the MP systems’ theory for modelling metabolic processes. In this paper, a new algebraic formulation of inverse dynamical problems, based on MP grammars and Kronecker product, is given, for further motivating the correctness of the LGSS (Log-gain Stoichiometric Stepwise) algorithm, introduced in 2010s for solving dynamical inverse problems in the MP framework. At the end of the paper, a section is included that introduces the problem of multicollinearity, which could arise during the execution of LGSS, and that defines an algorithm, based on a hierarchical clustering technique, that solves it in a suitable way. Key words: Metabolic P systems, dynamical systems, dynamical inverse problems, Kronecker product, stepwise regression. 1 Introduction Metabolic P (MP) systems are a particular class of cell-like P systems [33, 34, 36, 35] introduced by Vincenzo Manca in 2004, for modelling metabolic processes [29]. An MP system is essentially a particular type of deterministic discrete dynamical system which inherits from the P systems’ framework a native similitude with the functioning of a living cell. MP systems share with P systems the multiset rewriting mechanism as their fundament. However, while P systems are essentially unconventional computational models, MP systems are intended to generate dynamics instead of computations. Namely, their aim in modelling biological phenomena is that of finding the multiset rewriting mechanism underlying an observed biological behaviour. Metabolic P systems can be considered as the result of a research activity initiated in 1990s with some initial works [15, 28, 30]. They are different, with respect to other “P variants” applied in the context of systems biology [3, 4, 6, 38, 39]. The main difference is in their determinism. In fact, their basis are MP grammars, where multiset transformations are regulated by functions in a 2 V. Manca, L. Marchetti deterministic way [19]. An MP system is an MP grammar equipped with a temporal interval τ, a conventional mole size ν, and substances masses, which specify the time and population (discrete) granularities respectively [19]. An MP grammar Gcan be considered as a generator of time series, determined by the following structure (n, m ∈N, the set of natural numbers): G= (M, R, I, Φ) where: 1. M={x1, x2, . . . , xn}is a finite set of elements called metabolites, or substances. A metabolic state is given by a list of nvalues, each of which is associated to a metabolite. 2. R={αj→βj|j= 1, . . . , m}is a set of rules, or reactions, with αjand βj multisets over Mfor j= 1, . . . , m. 3. Iare initial values of metabolites, that is, a list x1[0], x2[0], . . . , xn[0] providing the metabolic state at step 0. 4. Φ={ϕ1, . . . , ϕm}is a list of functions, called regulators, one for each rule, such that, for 1 ≤j≤m, and for some kj(0 ≤kj≤n) ϕj:Rkj→R. An MP grammar Gis parametric, when a set Pof parameters is added to G, and metabolic states include also elements of P(to which, the state assigns real values), therefore regulators may include parameters as their arguments . If G is parametric, also the time series of parameters has to be provided in order to specify G. An MP grammar can be easily representable by an MP graph [22]. Moreover, the set of the rules of the system can be also represented by a stoichiometric matrix A, which gives a sort of “matrix-like representation” of the stoichiometry (see Figure 1). An MP grammar Gdefines, for any x∈M, a time series (x[i]|i∈N, i > 0) in the following way. Let s[i] = (x1[i], x2[i], . . . , xn[i]) the (row) state vector of Gat step i, which can be seen as a function from the set of metabolites to R, then the flux ϕj(sj[i]) of rule rjat step i, is given by applying the regulator ϕjto sj[i], a substate of s[i] associated to rj, and constituted by kj components called the tuners of rj. If we consider the rule r2of the MP grammar given in Figure 1, for example, then the flux at step iis calculated by: ϕ2(s2[i]) = ϕ2(A[i], B[i]) =c2A[i]2+c3B[i]3 Inverse Dynamical Problems: An Algebraic Formulation Via MP Grammars 3 r₄ A B C c₁A² c₂A²+c₃B³ c₇+c₈C c₄B c₅+c₆AC² MP grammar Rules Regulators r₁: ∅ → A φ₁ = c₁A² r₂: A → B φ₂ = c₂A²+c₃B³ r₃: B → ∅ φ₃ = c₄B r₄: A → C φ₄ = c₅+c₆AC² r₅: C → ∅ φ₅ = c₇+c₈C c₁, c₂, ... , c₈ ∈ ℝ A[0], B[0], C[0] ∈ ℝ r₃ r₁ r₅ r₂ MP graph 1 -1 0 -1 0 0 1 -1 0 0 0 0 0 1 -1 ( ) 𝔸 = Stoichiometric matrix Action of r₂ on substances A, B and C (it consumes A and produces B) Action of r₁, r₂, r₃, r₄ and r₅ on substance B (produced by r₂ and consumed by r₃) } } a column for each rule r₁ r₂ r₃ r₄ r₅ A B a row for each substance C # # # # # Fig. 1. An example of MP grammar (where ∅denotes an empty multiset and substance symbols occurring in regulators denote the corresponding substance quantities), the stoichiometric matrix Ais directly deduced by the MP grammar on the top left corner. The MP graph on the top right corner is obtained by translating the rules in the source-target-edge notation [26]. where c2, c3are given real constants, and Aand Bare said to be the tuners of the rule r2. The value of x[i+ 1], for each x∈M, is given by the following equation, where αj(x) and βj(x) denote the multiplicities of xin the multiset αjand βj, respectively: x[i+ 1] = x[i] + m X j=1 [(βj(x)−αj(x)) ·ϕj(sj[i])]. More generally, if we denote by Athe stoichiometric matrix of the system and by Φ[i] = (ϕ1(s1[i]), ϕ2(s2[i]), . . . , ϕm(sm[i])) the row vector of fluxes at step i, it can be proved that [16]: (s[i+ 1] −s[i])T=A×ΦT[i] (1) that is, by transposition: s[i+ 1] −s[i] = Φ[i]×AT.(2) 4 V. Manca, L. Marchetti These last two equations define equivalently the Equational Metabolic Algorithm (EMA). In the following, the MP dynamics we will present are computed in MATLAB1by applying EMA. We refer to [19, 20, 18, 21] for a comprehensive presentation of the MP theory. The dynamics which can be modelled by MP systems can be very complicated even by considering simple MP grammars (i.e. with few substances and linear regulators). In [23] MP systems were successfully applied to the field of real periodical function approximation. The complexity of the dynamics compared to the simplicity of the MP grammar which calculates it by EMA, suggests that MP system theory can be a suitable framework for modelling biological dynamics. The procedure introduced in [23] to define the models has been widely extended in [25, 26] for defining the LGSS (Log-Gain Stoichiometric Stepwise) algorithm, which derives MP grammars generating time series of observed dynamics. LGSS can be applied independently from any knowledge about reaction rate kinetics and it represents the most recent solution, in terms of MP systems, of the dynamical inverse problem, that is, of the identification of (discrete) mathematical models of an observed dynamics and satisfying all the constraints required by the specific knowledge about the modelled phenomenon. The LGSS algorithm combines and extends the log-gain principles developed in the MP system theory [16, 17] with the classical method of Stepwise Regression [7], which is a statistical regression technique based on Least Squares Approximation and statistical F-tests [5]. LGSS has been implemented by Luca Marchetti in 2010 as a set of MATLAB functions. We refer to [24, 31, 32, 27] for some successful applications of LGSS and MP systems for discovering the internal regulation logic of phenomena relevant in systems biology. The starting point of the LGSS algorithm was the search for the right regulators associated to the reactions of an MP grammar which provide the observed time series when dynamics is computed by means of EMA. If we consider the role of each regulator, we realize that it affects the variations of many substances. Therefore regulators are constrained to satisfy altogether, at each step, an algebraic system based on the stoichiometry of the observed phenomenon. The crucial point for regulator determination was a special kind of regression formulated as “stoichiometric expansion” of EMA by means of an initial set of basic functions called regressors. In the next section we will introduce a new algebraic formulation of the stoichiometric expansion, based on MP grammars and Kronecker product, which better describes and motivates its adoption in LGSS for solving inverse dynamical problems. 1See http://www.mathworks.it/index.html for details on the MATLAB software. Inverse Dynamical Problems: An Algebraic Formulation Via MP Grammars 5 2 Stoichiometric expansion Given a system with nvariables x1, x2, . . . , xn, let us suppose to know the time series of these variables along time points 0,1, . . . , t. Let s[i] = (x1[i], x2[i], . . . , xn[i]) the (row) state vector at time i, and xj[i+ 1] −xj[i] = ∆j[i] for j= 1,2, . . . , n, then s[i+ 1] −s[i] = (∆1[i], ∆2[i], . . . , ∆n[i]) whence, from equation (2), we get Φ[i]×AT= (∆1[i], ∆2[i], . . . , ∆n[i]).(3) For the determination of the regulators which provide the best approximate solution of the system (3), which has munknowns (the mcomponents of the flux vector Φ[i]), LGSS applies a procedure called stoichiometric expansion. Let us assume that the regulators we are searching for can be expressed as linear combinations of some basic regressors g1, g2, . . . , gdwhich usually include constants, powers, and products of substances, plus some basic functions which are considered suitable in the specific cases under investigation: ϕ1=c1,1g1+c1,2g2+. . . +c1,dgd ϕ2=c2,1g1+c2,2g2+. . . +c2,dgd(4) . . . =........................... ϕm=cm,1g1+cm,2g2+. . . +cm,dgd. Let us consider the t-expansion Φt 1, Φt 2, . . . , Φt mof regulators as the vectors constituted by the right members of equations (4) evaluated along tsteps (where the values of all the variables of the system are supposed to be known): Φt 1=c1,1Gt 1+c1,2Gt 2+. . . +c1,dGt d Φt 2=c2,1Gt 1+c2,2Gt 2+. . . +c2,dGt d(5) . . . =.............................. Φt m=cm,1Gt 1+cm,2Gt 2+. . . +cm,dGt d. Now, let Cd 1, Cd 2, . . . , Cd mbe the unknown column vectors, of dimension d, constituted by the coefficients of the regressors providing the linear combinations of regulators ϕ1, ϕ2, . . . , ϕmwe are searching for, and C= (Cd 1, Cd 2, . . . , Cd m) 6 V. Manca, L. Marchetti the matrix having these vectors as columns. Moreover, let ∆t 1, ∆t 2, . . . , ∆t nbe the column vectors of dimension tconstituted by substance variations of substances, from step ito step i+ 1, for 0 ≤i≤t−1, and ∆= (∆t 1, ∆t 2, . . . , ∆t n) the matrix having these vectors as columns. Let also Φtbe the following matrix constituted by mcolumn vectors of telements: Φt= (Φt 1, Φt 2, . . . , Φt m). Finally, let G= (Gt 1, Gt 2, . . . , Gt d) the matrix, of dimension t×d, having as columns the vectors obtained by evaluating the regressors g1, g2, . . . , gdon the tobserved time points. With the notation above, the system of equations (5) becomes: G×C=Φt.(6) Now, it easily follows from (3) that: Φt×AT=∆(7) where the exponent Tdenotes the matrix transposition. Therefore, by combining equations (6) and (7), we finally obtain the t-expansion of the system (3) as: G×C×AT=∆. (8) The coefficients of Care the unknowns which needs to be estimated by LGSS. We show now that they can be obtained by a Least Square Estimation deduced by equation (8), by using direct product ⊗between matrices, also called Kronecker product [9, 10, 40, 41], which results a special case of tensor product used in linear algebra and in mathematical physics. Given two real matrix A, B of dimension n×mand t×drespectively, then the direct product: A⊗B is the matrix, of dimension nt ×md, constituted by nm blocks Bi,j, such that, if A= (ai,j |1≤i≤n, 1≤j≤m), then Bi,j =ai,jB(in Bi,j all the elements of Bare multiplied by ai,j , see Figure 2). The Kronecker product is bilinear and associative, that is, it satisfies the following equations: A⊗(B+C)=(A⊗B)+(A⊗C) (A+B)⊗C= (A⊗B)+(A⊗C) (kA)⊗B=A⊗(kB) = k(A⊗B) (A⊗B)⊗C=A⊗(B⊗C). Inverse Dynamical Problems: An Algebraic Formulation Via MP Grammars 7 a b c d e f ⊗α β γ δ =       aα β γ δ bα β γ δ cα β γ δ  dα β γ δ eα β γ δ fα β γ δ        Fig. 2. An example of Kronecker product of two matrices. Moreover, matrix direct product verifies also the following equations: (A⊗B)×(C⊗D)=(A×C)⊗(B⊗D) (A⊗B)T=AT⊗BT (A⊗B)−1=A−1⊗B−1 where the exponent Tdenotes transposition and the last equation holds only when the involved matrices are invertible. Let us denote by vec(W) the vectorization of the matrix W, obtained by concatenating in a unique column vector all the columns of Win their order. Then, a general property of matrix direct product asserts that [10]: A×X×B=Y iff (BT⊗A)×vec(X) = vec(Y).(9) Therefore, if we apply equivalence (9) to equation (8) we obtain: (A⊗G)×vec(C) = vec(∆) (10) where the stoichiometric matrix is multiplied, by Kronecker product, with the regressor matrix and the result is multiplied with the vectorization of the regressor coefficient matrix, and then equated to the vectorization of the substance variation matrix, by providing nt equations with md unknown values. The system of equations given in (10) is the stoichiometric expanded system calculated by LGSS. According to the Least Square approximation method [43, 13], if nt ≥md, then the best approximation to vec(C), minimizing the difference between the two members of equation (10), is given by the following vector: (A⊗G)T×(A⊗G)−1×(A⊗G)T×vec(∆).(11) 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 xcannot exceed the quantity of x). In Figure 3 are represented the regressor matrix Gand the substance variation matrix ∆which are used by LGSS for least-squares approximating the coefficients c1, . . . , c8of the MP grammar given in Figure 1. 8 V. Manca, L. Marchetti Regressor matrix (A[0])² (A[1])² ... (A[t-1])² A² { B[0] B[1] ... B[t-1] { B C[0] C[1] ... C[t-1] { C A[0]·(C[0])² A[1]·(C[1])² ... A[t-1]·(C[t-1])² AC² { ) (B[0])³ (B[1])³ ... (B[t-1])³ B³ { [ ] t [ ] t [ ] t [ ] t [ ] t Substance variation matrix A[1] - A[0] A[2] - A[1] ... A[t] - A[t-1] ∆ ( { [ ] t A B[1] - B[0] B[2] - B[1] ... B[t] - B[t-1] ∆ { [ ] t B C[1] - C[0] C[2] - C[1] ... C[t] - C[t-1] ∆ { ) [ ] t C 1 1 ... 1 ( { 1 [ ] t Δ= 𝔾 = Fig. 3. The regressor matrix Gand the substance variation matrix ∆used for approximating the coefficients c1, c2,...,c8of the MP grammar given in Figure 1. However, the approximation given by (11) cannot in general be considered the best way for solving the inverse dynamical problem. In fact, apart the computational cost of considering all the dregressors at same time, several reasons suggest to follow a gradual strategy in the determination of a subset of regressors and their corresponding coefficient which provide the best approximation to the given dynamics. There are two main requirements which are essential for an appropriate application of least squares method: the linear independence among the regressor expansions and the parsimony of the set of regressors. In other words, the best approximation is obtained by determining a parsimonious set of linearly independent regressors ensuring an error under a given threshold. Linear independence is a requirement of least squares method and is solved by considering systems of equations which have been stoichiometric expanded. The parsimony of the model, instead, avoids problems of overfitting. In fact, the more regressors are considered in the model, the less is the degree of freedom left for the error [1]. This implies that solution fits very well with the dynamics on the observation points, but it is too constrained to them for behaving in a satisfactory way outside them (i.e. the model fits well the data, but it has not predictive power, see Figure 4 for an example). In order to cope with the requirements explained above, LGSS integrates the least squares approximation of stoichiometric expanded systems with a regression strategy based on a step-wise approach as defined in [26]. Such kind of approach permits to define the model, step by step, by inserting into the model only those expanded regressors (among the columns of the matrix given by the direct product A⊗G) which satisfy specific statistical tests. In this way, we can obtain MP models which fit the dynamics and that comprehends a small set of regressors. 3 Problems related to the regression in LGSS The stepwise approach adopted in LGSS is based on the assumptions which are at the basis of the classical multiple regression model [1]. These assumptions concern with some properties of the expanded regressors (i.e. they must be linearly Inverse Dynamical Problems: An Algebraic Formulation Via MP Grammars 9 independent and, possibly, not correlated2) and with the probability distribution of the errors associated to observations in considered time series (i.e. the errors should be normally distributed with mean zero). When one or more of these assumptions are not completely satisfied, some mistakes can occur in the definition of the regulators. In particular, there are several problems which we need to be aware of in the context of multiple regression. Some of them have been discussed in [26] and can be solved by substituting the ordinary least squares with other estimation methods based on the weighted least squares [42] or on the generalized least squares [12]. Here we focus on solving the problem of multicollinearity, which consists in having regressors that are highly correlated among them. This is the most common problem occurring in LGSS and also one of the most difficult to be solved [1]. When we develop a new MP model, we hope to have a strong correlation between each expanded regressor and the dependent variable vec(∆), but we do not want to have expanded regressors correlated among them. In fact, this phenomenon may cause errors in the selection of the right set of regressors during the execution of the stepwise regression. In the case of perfect collinearity, the regression algorithm breaks down completely (because the matrix given by the direct product A⊗G has not maximum rank). Since in LGSS usually regulators are assumed to be linear combinations of polynomial regressors, then it is very common to meet multicollinearity problems. 2The correlation between regressors is intended to be calculated by means of the Pearson’s correlation coefficient [37], which ranges from −1 to 1 and provides a measure of dependence between the behaviours of two magnitudes (−1: perfect anti-correlation; 0: no correlation; 1: perfect correlation). Fig. 4. Comparison between the predictive power of two regression models: a 13-degree polynomial b Y=c0+c1X+c2X2+. . . +c13X13 (depicted by the continuous line) and a least squares line (depicted by the dotted line). The dataset used to calculate the models are the 14 points depicted as blue circles, the last point represented by the red star is the value of Ywe want to predict with our models. The 13-degree polynomial is a perfect example of model which overfits the data: in fact, it provides a perfect fit for all the points of the dataset, but it completely fails the prediction of Yin the 15th data point. 16 V. Manca, L. Marchetti [36] Gh. P˘aun, G. Rozenberg, and A. Salomaa, editors. Handbook of Membrane Computing. Oxford University Press, 2010. [37] K. Pearson. Notes on the History of Correlation. Biometrika, 13(1):25–45, 1920. [38] F.J. Romero-Campero and M.J. P´erez-Jim´enez. Modelling gene expression control using P systems: The Lac Operon, a case study. Biosystems, 91(3): 438–457, 2008. [39] A. Spicher, O. Michel, M. Cieslak, J.L. Giavitto, and P. Prusinkiewicz. Stochastic P systems and the simulation of biochemical processes with dynamic compartments. Biosystems, 91(3):458–472, 2008. [40] W.H. Steeb. Matrix Calculus and Kronecker Product with Applications and C++ Programs. World Scientific Publishing, 1997. [41] W.H. Steeb. Problems and Solutions in Introductory and Advanced Matrix Calculus. World Scientific Publishing, 2006. [42] T. Strutz. Data Fitting and Uncertainty. A practical introduction to weighted least squares and beyond. Vieweg+Teubner, 2010. [43] J. Wolberg. Data Analysis Using the Method of Least Squares: Extracting the Most Information from Experiments. Springer, 2005.