Search :   Advanced Search 

Home   |   Join   |   Contact     

   
Journal Details
 
Article usage
Total views: 1140
[From(publication date):
December-2011- May 22, 2012]
Breakdown by view type
HTML page views : 329
PDF downloads : 441
XML downloads : 370
 
 
Subscribe Here
Enter your name :*
Enter your Email : *
 
 
 
 
Research Article Open Access
Towards an Exact Reconstruction of a Time-Invariant Model from Time Series Data
Michael A. Idowu1* and James Bown2
1Scottish Informatics Mathematics Biology and Statistics (SIMBIOS) Centre, School of Contemporary Sciences, University of Abertay, Dundee, Scotland, United Kingdom
2Institute of Arts, Media and Computer Games, University of Abertay Dundee, Scotland, United Kingdom
*Corresponding author: Dr. Michael A. Idowu
Scottish Informatics Mathematics Biology and Statistics (SIMBIOS) Centre, School of Contemporary Sciences
University of Abertay, Dundee
Scotland, United Kingdom DD1 1HG
E-mail: m.idowu@abertay.ac.uk
 
Received September 08, 2011; Accepted November 10, 2011; Published November 23, 2011
 
Citation: Idowu MA, Bown J (2011) Towards an Exact Reconstruction of a Time- Invariant Model from Time Series Data. J Comput Sci Syst Biol 4: 055-070. doi:10.4172/jcsb.1000077
 
Copyright: © 2011 Idowu MA, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License,which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
 
Abstract
 
Dynamic processes in biological systems may be profiled by measuring system properties over time. One way of representing such time series data is through weighted interaction networks, where the nodes in the network represent the measurables and the weighted edges represent interactions between any pair of nodes. Construction of these network models from time series data may involve seeking a robust data-consistent and time-invariant model to approximate and describe system dynamics. Many problems in mathematics, systems biology and physics can be recast into this form and may require finding the most consistent solution to a set of first order differential equations. This is especially challenging in cases where the number of data points is less than or equal to the number of measurables. We present a novel computational method for network reconstruction with limited time series data. To test our method, we use artificial time series data generated from known network models. We then attempt to reconstruct the original network from the time series data alone. We find good agreement between the original and predicted networks.
 
Keywords
 
Mathematical modelling; Dynamic systems; Reverse engineering; Inverse problems; System identification and parameter estimation; Time series analysis and Time-invariant models; Network inference algorithm
 
Introduction
 
Interaction networks may be used to describe a wide range of biological phenomena, spanning genetic [1] and biochemical networks [2], analysis of clinical orbiomolecular data [3], and personalised medicine [4]. These network-based approaches seek to describe dynamic systems as interaction networks, where the nodes in the network represent the measurables and the edges represent interactions between any pair of nodes. Additionally, networks may have weighted edges where the weights relate to the strength of the interactions [1]. For example, gene networks are commonly represented by directed graphs where the nodes of the graph are genes and the directed edges are causal relationships between genes [1]. In metabolic pathway modelling, the first step involves the development of appropriate functions that describe the behaviour of all the constituents of the system and identification of all components of interest with their symbolic names with (directional) arrows showing which components modulates the ows into, between, and out of components [2]. In such cases, the set of nodes is defined by the measurables; importantly, the edges and their weights, and so the network topology, must be inferred from the data.
 
Many systems, within and beyond biology, may be described by interaction networks. For example, within biology gene regulatory networks and intra-cellular signalling networks are essentially interaction networks: nodes represent genes and bio-molecular species respectively; edges represent interactions, which may change under different perturbations. Such biological networks are characterised by complex interactions and feedback loops [5]. Beyond biology, our socio-technical infrastructures are readily described by interaction networks, both in terms of their individual characteristics and their inter-network interactions [6]. Our energy, transport, water and telecommunications networks are a complex adaptive system [7] akin to biological interaction networks. These complexities, in biological and engineered systems, make system dynamics challenging to interpret. As a result, system-scale methods that support network identification and parameter estimation methods are important components in developing our understanding of these systems [8]. Simple, effective, reasonable and robust algorithms for constructing underlying network topologies are in demand for understanding the dynamics of complex systems and processes.
 
Time-evolution of complex processes may be represented in the form of time series data. These time series data are state vectors measured over time in linear order and often at regular intervals. Describing the dynamic properties of the observed time series of a system can help identify and model patterns, predict behaviour, and approximate the underlying process generating the data [3]. Systems that generate time series data may be highly complex and incapable of exact description, so approximating the underlying process through parametric models can often help capture key features of data relevant for the purpose of interest [9].
 
Here, we propose and present a novel computational, robust method (solution) for constructing weighted interaction networks based on time series data.
 
The construction of an interaction network from time series data may be recast as seeking to identify a mathematical model that relates a given time point to its successor, consistent with every time step and for all measurables. This model may be expressed in the form of a n by m matrix, for n time points and m measurables. For a given time point t, the mathematical model must relate a measured value xi at t, xi (t), to its value in the subsequent time point xi (t+1). Moreover, this relation is required for all x(1…m) and for all t(0… (n - 1)). This mathematical model, referred to here as a transformation matrix, must be calculated from the data, i.e. the model and model parameters must be inferred from data. This process may be described as an inverse problem, where we must identify the transformation matrix solution to the system describing the time series data.
 
Finding a time-invariant matrix to the time series inverse problem may require that the number of time points provided be at least equal to the number of measured variables +1. Whenever the available data set is limited in size, the number of time points is less than or equal to the number of measured variables, solution identification becomes more challenging.
 
The method proposed here requires a minimum of 3 time points to be able to reconstruct a network. As we show later, as the number of time points increases the network reconstruction process becomes more accurate. However, even with (only) 3 time points, we show that the method produces a data-consistent model of a given data set. A model is data-consistent if it is capable of reproducing exactly the original time series data that was used to infer the model.
 
Any time series data set can be described as a system of ODEs, whose variables are the measurables in the physical system. For example, a gene regulatory network with multiple genes can be represented by a set of nonlinear ordinary differential equations (ODE) with the expression level of each gene as a variable [10]. Reverse engineering of such networks, through gene expression data analysis and reconstruction of gene regulatory networks, involves revealing the underlying network of gene-gene interactions from the measured dataset of gene expression. This process can involve some form of mapping of the observed data to a reconstructed and representative network model inferred from data using reverse engineering techniques [3]. With respect to the identification of the transformation matrix, we show below how variables in the ODEs model map to matrix elements and our algorithm operates on this matrix to reverse engineer interaction network topology and edge weightings. Solving this inverse problem has implications for data modelling in systems biology generally and in particular in personalised medicine through for example protein interaction network modelling.
 
Our algorithm searches for a solution to the inverse problem, and the algorithm identifies the same solution every time for a given problem. However, if the size of the available dataset is small, i.e., total number of time points < number of measured variables, other potential solutions may still exist. Our approach is different from other forms of optimisation algorithms such as genetic algorithms that are based on finding the best set of model parameters within a search domain, because starting points for the domain search process are decided prior to parameter estimation. Here, no single parameter is fixed prior to estimation and no start point comes into play; the solution is purely based on experimental data. Though some parameters can be set to fixed values, our parameter estimation technique is purely based on and driven by experimental data, without any need for fixed parameters.
 
We outline how dynamic systems may be described by systems of ordinary differential equations, and present a matrix-based approximation to the solution. We then describe our algorithm to solve the inverse problem through identification of the values of elements in this matrix. We demonstrate the capability of the algorithm to reverse engineer interaction networks with a worked example, provide summary statistics on algorithm performance and comment on the uniqueness of the solutions discovered.
 
Ordinary Differential Equation Systems
 
It is common to use ODE formalisms to model and analyse data in complex dynamical systems [10]. Sometimes, simple linear ODE models are sufficient representations of the system in order to identify essential relationships among network components based on time series data analysis. Solutions to systems of first-order ODE may involve matrix factorisation to approximate a matrix exponential [11].
 
Specifically, the solution of the homogenous equation
 
                                                                                                                         (2.1)
 
where X is a column vector representing the dependent components of the system, t is a time variable and J, the relative rate of change of X(t), is a matrix often referred to as jacobian, is given by
 
                                                                                                                           (2.2)
 
Where X(0) = X0 the initial state of the system for a system of linear differential equations, involving a transforming matrix J and the solution [ej*t * X ] a function of (J, X, t), where X is a column vector representing the dependent components of the system, and t is a time variable that represents a regular time interval between any two successive states, if J (the transformation matrix inferred from time series data) is derived such that it remains unchanged throughout the system (from the initial condition to steady-state, and after), then the model X = J*X may be said to be time-invariant because J is not dependent on time. In other words, the parameters of J do not change and are not functions of time.
 
Let Xt and Xt+1 be the state vectors known from a given time series data, and assume that the time interval tc is regular and of small magnitude. Then , and the inverse problem is mathematically equivalent to
 
                                                                                                                (2.3)
 
There may be more that one J matrix that satisfies the above equation, so the primary objective is to find the best J that describes - with least error - the transformation from any state Xt to the next state Xt+1.
 
Relative rate of change
 
J, the relative rate of change in X(t), may be described as
 
                                                                                                                                  (2.4)
 
that is, J is the absolute rate of change X in relation to the present state value X. Since it is true that therefore one may describe J in terms of which means that J itself is that relative rate of change, and it is equivalent to the rate of change in the logarithm of X(t) [12].
 
The jacobian, in a sense, can be thought as a representation of the rate and extent of change over time for any component and its temporal influence on other components. The jacobian matrix may thus be regarded as a construct for describing system dynamics. The mathematical significance of determined eigenvalues, or chacteristic values, of the jacobian can help inform understanding of the behavior of systems near a stationary point. By observing the eigenvalues of the jacobian matrix in a given neighbourhood, an indication of the stability of the system can be obtained. If the eigenvalues all have a negative real part, the system is said to be stable [13]. A positive real part in any of the eigenvalues means the system may be unstable at that point, exhibiting large transient peaks [13]. Regarding the determinant of the jacobian matrix, the absolute value of the determinant of the (jacobian) transformation matrix indicates the overall rate of change of all the measured variables [14]. A necessary and sufficient condition for the jacobian matrix invertibility is that the magnitude of its determinant > 0 [14].
 
The algebraic representation of an ODE-based model solution for a time series problem is simple and straightforward.
 
                                                              (2.5)
 
And this is related to the eigenvectors and eigenvalues representation as follows
 
 
where the parameters λi and υi represent the eigenvalues and eigenvectors, respectively; and the initial condition, [X1(0),X2(0),…Xn(0)] is favourably chosen, i.e., decomposed and approximated, as the linear combination of the eigenvectors of the jacobian matrix using the parameter set [p1,p2,… ..,pm]. Each measurable Xi (component within a network) is a function based on the initial condition. The initial condition (first measurement) or any state vectors at any timepoint may be rewritten (decomposed) in a form such that the parameter set [p1,p2,…..,pm] becomes fixed [12].
 
Complex nonlinearity in systems of nonlinear ODEs
 
Since most complex systems are nonlinear in nature, nonlinear models are required to describe them fully. Moreover, some systems may require that second-order ODEs be used to formulate a sufficient description of their dynamics due to the complex nature of the processes involved. However, in such cases, the inference of model parameters becomes a much more difficult task. In particular, the appropriate model structure needs to be identified, or at least estimated with a reasonable degree of confidence, before model parameters are estimated.
 
Here, we suggest that first-order linear ODE solutions based soley on time series data, i.e., with no assumptions made about network structure, may be sufficient to derive valuable information about the most important links among variables if good search algorithms for network inference are employed.
 
In some cases modellers predetermine network structure prior to parameter estimation, e.g. [15]. By keeping model structure fixed in this way, the inference procedure is restricted and identification of other solutions that do not adhere to the fixed structure constraint is not achieved. Such techniques themselves are thus limited since they are only successful when the constraint imposed on the system is met. Of course, the rationale for fixing the network topology is that it eases identification of potential solutions for underdetermined systems [15]. As a result, many modellers tend to assign some fixed values to subsets of model parameters to ease the fitting process. A good algorithm is one which can accommodate fixed parameters and at the same time does not require these to find an accurate solution.
 
We propose such an algorithm, effected via a combination of pre-conditioning, regression, analytical techniques, half and or S-system approaches based on time series data. Our rationale for proposing linear ODE solutions to solving inverse problems in time series are as follows:
 
(1) Second-order differential equations can be recast into first-order (ordinary) differential equations;
 
(2) Any nonlinear ODE can be cast or approximated into a power-law form called (S-systems or half systems);
 
(3) Nonlinear half systems are equivalent to systems of log-linear differential equations;
 
(4) Log-linear differential equations may be used using the same techniques for solving systems of linear differential equations.
 
Therefore, the development of a robust method for solving systems of linear differential equations is relevant to non-linear problems. Multiple regression and reverse engineering techniques, including the logarithmic inverse function, are important steps we have considered. Here, however, we limit the focus of the paper to time series analysis under systems of linear ordinary differential equations.
 
Methodology
 
In a system of linear differential equations an inverse problem may be written as in (2.4):
 
where and X are known vectors of same length, and J is the unknown matrix that must be identified. Note that there is difference between a general system of n linear differential equations with unknown (jacobian) matrix parameters and a general system of n linear equations with unknown vector parameters. The latter is much simpler to solve due to the reduction in the number of unknown parameters. However, the formulation of the inverse problem remains the same in structure.
 
A general system of m linear equations with m unknown parameters is of the form b = A.x (3.1)
 
with the following matrix equations:
 
                                                                            (3.2 )
 
                                                                                     (3.3 )
 
where b=[b1,b2,…,bm]T are the constant terms, a11,a12,…,amm are the coefficients of the system, and x=[x1,x2,….,xm]T are the unknown parameters; note [υ]T denotes the transposed vector [υ]. Finding a solution to the system above involves searching for values in the x vector where all the m equations are simultaneously satisfied and valid.
 
A general system of m linear differential equations with an unknown m x m jacobian matrix J is of the form and has the following matrix equation:
 
     (3.4) 
 
where = X(t) is a known state vector recognized as the tth vector of X (not X at time t), X is the derivative vector which may be calculated from two known state vectors (Xt and Xt+1) as where tc is the interval of separation, and J is the unknown m x m jacobian matrix. Therefore, at least two state vectors of the same length are required to define an inverse problem in systems of linear differential equations. The following multi-state representation defines an inverse problem involving n+1 different states:
 
                 (3.5)
 
which is equivalent to
 
          (3.6)
 
assuming the states values are captured at regular time interval tc. The smaller the value of tc the better the outcome of this derivative approximation; this is because the linear approximation, of improves as tc gets smaller and closer to 0 [12]. The solution to this system of linear differential equations:
 
 (3.7)
 
Is thus
 
                   (3.8)
 
Consequently solving an inverse problem in a system of linear differential equations is equivalent to identifying that fits the data best. This requires optimal estimation of the parameters of J; calculating all the entries (parameters) of the matrix of all first-order partial derivatives of vector-valued functions of variable i with respect to another variable j, where Xi or Xj represents the variable function of component i or j respectively. Using simple variables, the solution may be rewritten as: which is on the one hand a representation of system of linear equations, X(N+1)=E.X(N), and on the other hand a direct interpretation of a system of nonlinear equations, , implying that . Consequently E is equivalent to the matrix exponential (function) of the matrix product J.tc. The time constant tc is easily calculated as the difference in time between T1 at any state in X(N+1) and T0 its previous state in XN. E can easily be approximated by our new regression techniques (described below), and these are variant forms of regression for ODE systems. Often regression analysis is used to understand the relation between two or more interrelated variables. In systems of linear differential equations, regression analysis can be used to infer causal relationships or transformations between the model variables and states. So in principle, the state transformation (or the transposive regression) matrix, E, is the matrix exponential of (J.tc).
 
From these definitions, different approaches in analysis may be followed. One which is of worthy mention is the logarithm of E, which in this case is the actual result (J.tc) being derived from E such that its exponential equals E.
 
Approximating a value of E may comprise using a regression technique having the steps:
 
1. Acquire time series data with the number of time points ≥3;
 
2. Preprocess the measured state values of one or more components of the system using matrix transposition;
 
3. Undertake regression analysis of the resultant data using a Moore-Penrose pseudoinverse technique;
 
4. Postprocess the resultant data using matrix retransposition;
 
5. Calculate the logarithmic inverse of the retransposed result through application of eigenvectors and eigenvalues techniques;
 
6. Scale down the resultant data by factor (magnitude) of the time interval used.
 
The implementation of the steps 2-6 is shown in Sections 3.1.1 and 3.1.2 below. Section 4 considers the provision of time series data (step 1).
 
Transposive and repressive regression methods
 
In solving a system of linear differential equations, we define the solution to an inverse problem as one conditioned on the property:
 
 
Or simply:
 
 
where , the time interval tc is assumed to be regular; J is unknown at this point and must be identified, and X(t) and X(t+1) are known arrays of column vectors, each a representation of system states at two successive time points, termed before and after. In this paper we present for the first time two new algorithms to derive J, namely:
 
1. Transposive Regression Algorithm
 
2. Repressive Regression Algorithm
 
Derivation of the transposive regression algorithm:
 
Steps
 
(1)
 
(2)
 
(3)
 
(4)
 
(5)
 
(6)
 
In steps 1-2 recasting the problem by matrix transposition is essential, because each state is represented by a column-vector either in X(before) or X(after), where X (before) is an array of states before the transformation X (before) = [X (0) X (1)… X (t-1)] and X (after) is an array of states after the transformation X (after)= [X (1) X (2)… X (t)]
 
Steps 3-4 illustrate an application of the Moore-Penrose pseudoinverse, a widely known type of matrix pseudoinverse, independently introduced by Moore [16], Bjerhammar [17], and Penrose [18]. Finally, in steps 5-6, retranspositions put E1 in proper order.
 
Derivation of the repressive regression algorithm:
 
Steps
 
(1)
 
(2)
 
(3)
 
(4)
 
(5)
 
(6)
 
(7)
 
(8)
 
Here, in step 1 we first repress the equation by subtracting X(before) from both sides of the equation. In steps 2-3 we recast the problem by matrix transposition. In steps 4-5 the Moore-Penrose pseudoinverse is applied. And the re-transposition step is introduced in steps 6-7. The identity matrix (I) is added to both sides of the equation to derive E2 on the left hand side in step 8.
 
The search for the jacobian matrix solution
 
It is not difficult to calculate the jacobian matrix once either E1 or E2 is found. There are two different methods to consider:
 
(1) using eigenvalues and eigenvectors;
 
(2) using a new approximation technique, presented here for the first time.
 
The difficulty in finding J is in calculating the principal matrix logarithm of E1 or E2; that is, the exact inverse of
 
 
 
Where logm (. . . ) represents the matrix logarithm function.
 
Application of eigenvalues and eigenvectors
 
Assuming that E1 or E2 is diagonalisable, the following method may be used to obtain the jacobian matrix from E1 or E2. First we seek to find the matrix υ of eigenvectors of E1 or E2 as appropriate, referred to here as Em for convenience (for either case). Each column of υ is an eigenvector of Em. We then find eigenvectors of Em from υ and Em as . Next we replace each diagonal element of eigm by its natural logarithm and calculate the natural logarithm of Em as . Finally, J is then calculated from logm(Em) as follows:
 
                                                                    (3.9)
 
Note, only real values for parameters of J are considered.
 
A new approximation technique for calculating the matrix logarithm
 
It is generally known that is approximately equal to on condition that is a number and small enough [12]. Here, we will introduce a scaling factor μ to tc in order to approximate such that . We may calculate the value of J to be:
 
                                                                                                   (3.10)
 
                                                                                                   (3.11)
 
                                                                                                              (3.12)
 
                                                                                                                      (3.13)
 
By substituting a small value for μ in the above equation, e.g. 10-4≤μ≤10-11, we may approximate J. Note, the smaller the magnitude of this value the better the result of approximating J becomes. However, care should be taken not to allow the magnitude of this value to be smaller than this range in order to ensure that it is not approximated to zero internally. Note this new approximation technique is equivalent, in terms of the solution obtained, to (3.9).
 
Linking jacobian matrices and network models
 
The inference of model parameters from experimental data sometimes requires multivariate regression to be performed on experimental data with the aim of minimising the residual error between the model and the data, particularly in systems of linear and nonlinear ordinary differential equations (ODE). For example, using the following system of linear ordinary differential equations:
 
 
one can interprete the jacobian matrix to mean a direct representation of the network interaction and systems dynamics as indicated in (Figure 1).
 
Figure 1: Structure representation: (a) jacobian matrix of a network; (b) its corresponding network with inter-connected nodes.
 
Results
 
Method Validation
 
To test our approach for reconstructing network models from time series data, we use artificial data generated from known network models. This data is generated by simulating time series data from those models, and reconstructing the original network from the time series data alone with no knowledge of the generative model. Importantly, the original network is not provided to the reconstruction method - only the time series data, and this provides a source of independent test data not used in method construction. We test our method on 100 randomly generated networks as explained in section 4.5. We consider two data discretisation approaches for generating time series data, namely:
 
(1) discretisation using simple continuous models;
 
(2) discretisation using application of eigenvectors and eigenvalues.
 
Our objectives are:
 
(1) to discretise (multivariate time series) data with known network models (jacobian matrices);
 
(2) to ensure that the simulated data is noise-free (noiseless);
 
(3) to ensure that important features of the data such as correlation between the variables are preserved;
 
(4) to test and evaluate the performance of our inference algorithms based on minimum number of states ( <= number of measured variables +1).
 
The test methods avoid independent simulation of time series data of any single variable; only multivariate discretisation is used. To demonstrate the performance of our inference algorithms we also avoid using the integral function during the discretisation process. First we describe the discretisation techniques used to generate data. We present an expanded example of the method for each data type. Finally, we present summary statistics of performance of the method for reconstruction of a large number of networks generated at random.
 
Data discretisation using a simple continuous model
 
The solution to a system of ordinary differential equation is X(t) = eJtX0 as noted previously, so in order to discretise at regular time step intervals t (in this example our time step is 0.25 seconds), we define our discretisation process at any state k as :
 
                                                                                       (4.1)
 
Or at state k+1
 
                                                                                         (4.2)
 
Where X(0), X(k-1), X(k), X(k+1) are the vector values at the time points 0, k-1, k, and k+1, respectively. 4.2.1. Example 1. Define J as
 
 
and the initial condition, X(0), as:
 
 
Therefore, the state vector is:
 
 
where t=0.25 seconds. Likewise, state vector is then calculated and the result is:
 
 
Note that X(6) could have been calculated from X1 as: . If we defined a time series dataset DS1=[X(1) X(2) X(3)X(4) X(5) X(6) ], i.e., a time series data set for the next six states at regular interval of 0.25 seconds after the initial condition, then DS would be:
 
 
Data discretisation using application of eigenvectors and eigenvalues
 
It is widely known that the solution set of the system of ordinary differential equations can be represented by any combination of exponential functions of eigenvalues and their eigenvectors in the form:
 
 
where the initial condition, X0 , is represented as a linear combination of the eigenvectors of J [12]. Through further analysis, we introduce the matrix form of the initial condition as
 
 
so that we might present our general matrix form solution to be
 
 

where each column vector in
 
 
is an eigenvector of the jacobian matrix J and the parameter set {λ1, λ2,.λn} is a set of the eigenvalues of J. The parameter e is used to denote the exponential of 1. Note that the parameter set {p1, p2,.pm} can easily be calculated by regression analysis [12]. Using the example in the previous section where
 
the eigenvalues and eigenvectors of J are calculated to be:
 
real(eigVec) =
 
 
Imag(eigVec) =
 
 
eigVec=
 
 

eigVal=
 
 
Setting the initial condition to
 
 
implies that
 
 
with this estimated parameter set found through regression analysis. Therefore, we define our second discretisation process at any time point t as :
 
 

Therefore, the timepoint at time t 0.25 secs becomes:
 
 
We define time series dataset DS2=[X(1) X(2) X(3)X(4) X(5) X(6) ]
 
DS2 =
 
Reverse engineering results
 
We now demonstrate how well our reverse engineering (inverse problem solution) algorithms work using the time series created in the last section where
 
 

and  we define

 
 
 
and tc = 0.25
 
4.4.1. Example 1a: Application of the Transposive Regression Algorithm. First calculate E1 from
 

 
Then use either of the matrix logarithm techniques, introduced in Section 3.4, to reverse engineer J
 
 
4.4.2. Example 1b: Application of the Repressive Regression Algorithm. E2 is calculated from
 

 

Reverse engineering J from E2 then produces:
 
 
Performance evaluation of algorithms
 
Network structure identification: The performance of the two algorithms introduced in this paper are evaluated in terms of their ability to identify original (unseen) network structures. Before parameter estimation, model structures should be determined. The approximated network structure is easily derivable from the jacobian matrix (derivable from the zero and non-zero entries in the matrix representation of the system of ODEs), regardless of the magnitude of the parameter entries. Therefore, we simplify a weighted jacobian maxtrix into its Boolean form, revealing the network structure (model architecture) in terms of presence and absence of links. These structures are in form of matrices containing only Boolean (0 or 1) entries: an entry of 1 depicts presence of an association (non-zero parameter in the jacobian matrix); otherwise 0 for no association. This initial approximation determines the network topology (see Figure 2).
 
Figure 2: Relation between the derived Boolean representation of the jacobian matrix and the corresponding network topology with inter-connected nodes.
 
Network connectivity ratio: Given a (Boolean) matrix representation of the inferred jacobian matrix (network structure), we define network connectivity ratio as where NumOfOnes and NumOfZeroes indicate the total number of 1s and 0s in the Boolean matrix respectively.
 
Metric criteria: Three key properties of network representation are used for performance evaluation:
 
a) ratio - a relative ratio of the number of missing correct links in the inferred matrix to the total number of links inferred; b) relative ratio in terms of number of incorrect links in the inferred matrix to the total number of links inferred; and c) the norm score based on and ratios - an indication of the degree of closeness of the predicted network to the original network, measured as
 
Algorithm performance
 
We tested our algorithms on 700 simulated datasets with a range of numbers of time points (4-10), all generated from networks with 25% connectivity ratio.
 
Well-formed jacobian matrices of artificial systems are required to generate data for performance evaluation of our method. A method for constructing nonsingular matrices was used to generate our artificial network data. Based on pameterisation of matrices implied determinants and minors, we were able to randomly generate a large set of new nonsingular jacobian matrices that were used to simulate test data by operating and manipulating products of factors of nonsingular matrices to guarantee that the jacobian matrices produced were not defective. With hundreds of such matrices, we were able to generate a wide range of different artificial data to test our method. The number of non-zero parameters in each of those matrices determined the network connectivity ratio for that system. Here, we fixed this to be 25% of the total parameters in a given matrix. Hence our matrices have 25% connectivity. Finally the results of the predicted networks were then compared to the corresponding information of the original networks recorded in the database.
 
Based on the performance evaluation criteria , and the rank (norm) score, we analysed the results of the inferred network structure and measured deviation of network structure size in terms of network connectivity percentage. We assume that an inferred network structure has the potential to be a reasonable or good result if its connectivity ratio is between 20 and 30%. Table 1 shows that at least 60% of the number of dependent variables (here 6) is the minimum number of time points required to obtain a reasonable or good inference result. However, with our methods the reconstructed models are often data-consistent, that is, have the capacity to simulate or reproduce exactly the given dataset, irrespective of the number of time points (even when there are only a few time points). The evaluation criteria have been established to measure variation in performance depending on the number of time points. The results confirm that performance improves with an increase in the number of time points (a decrease in rank score indicates an increase in performance). The whole performance evaluation process is automated and does not need supervision or user intervention.
 
The network connectivity ratio indicates the estimated number of model parameters identified from (and used to explain) the available data set, which means that the richer the dataset, in terms of number of time points, the better the probability of identifying the correct links (or model parameters) that are suitable or valid, e.g. in (Table 1), the first row indicated that on average at least 13 model parameters (predicted from data to be non-zero) were ascertained to be valid parameters out of a total number of 25, whereas the last row showed that at least 24 model parameters were identified as being valid. All (10) diagonal elements were included as valid entries by default.
 
Table 1: Summary statistics of algorithm performances.
 
Not every identified link is valid, although often the majority of them are. The network connectivity ratio reveals the complexity of the predicted networks, i.e.the total number of correct (and incorrect) links in the predicted network. The value, as previously noted, reflects the number of correct links predicted. As mentioned previously, the norm score based is based on the calculation . We assume that the lower the norm score the better the performance of the inference method. A high valued Norm Score might still produce a data-constistent model, but the jacobian matrix of such system would be too sparse for the output structure to be a reasonable representation of the underlying data. The result shows that performance, in terms of network structure identification, improves with an increase in the number of time points. The rank (norm) score (error ratio) value, a function of and ratios, confirms approximately success rate in inference of network structure for datasets with size of 10 time points. It is remarkable that with such limited information on data and no information on topology, our inference methods are able to infer completely (100%) the actual network structure at times.
 
In summary, the results show an improvement in method performance with increasing data size. Overall performance is close to optimum, i.e the original network is recovered with approximately 88.5% success rate on average, when the number of time points available is equal to number of measured variables even though this network is unseen to the algorithm and there are many possible data-consistent weighted networks. The main challenge is in keeping both the and ratios as low as possible. These two metrics may also be used in robustnessness and sensitivity analysis of any proposed method, keeping in mind that the primary objective of any proposed method is to minimise and values. Ultimately, the challenge is to preserve data-consistent model generation while at the same time maximising the likelihood of identifying the original set of links by inference.
 
Conclusion
 
Clearly, network structure identification and parameter estimation of dynamical systems are necessary steps in representing system dynamics in terms interaction networks. We demonstrate that algorithmic analysis of time series data may produce data-consistent models. On a promising note, the novel inference algorithms presented in this paper are identified, through simulation experiment and testing, to develop such dataconsistent models. As demonstrated in this paper using a worked example, there is a strong theoretical basis for their use in time series data analysis. Moreover, their utility is demonstrated by their performance result under testing conditions using artificial data sets generated in silico.
 
We assessed the performance of our unsupervised inference algorithm using hundreds of test networks, that is networks that were used to generate randomly valued, independent test data, through two different methods and similar in form (but not values) to the worked example, and importantly the underlying network structure was never provided to the algorithm in this validation. We demonstrated significant improvement in network reconstruction as more data became available, here increasing time series time points from 3 to 10, and showed very good performance as the data size tends to 10 time points. We recognise that 3 data points is a very small data set, but show that even with this limited time series data we are able to reconstruct a data-consistent network. Our algorithms are aimed at simplifying and standardising the methods of finding unique solutions whenever they exist and using those standardised methods to adequately find other potential data-consistent solutions in non-unique scenarios. Of course, as time series data reduces in size, the number of possible networks able to explain the data increases. Note that this ability to work with limited data can be combined with the capacity for the approach to include a priori knowledge, and this knowledge may substantially reduce the solution space. Consequently the approach can blend available knowledge with knowledge gaps to produce data-consistent models of system dynamics.
 
The application of this algorithm in systems biology extends to any time series data. Where those systems are better described by log-linear processes, our algorithmic approach is still valid. We are especially interested in the impact of cancer drug intervention strategies on the (human) cell signalling network (for example [20]). This cell signalling model describes the PI3K/AKT signalling network and considers the effects of different perturbations on the network response to growth receptor inhibition. By perturbing the system with various mutations, distinct regimes of functioning were observed in the network. Specifically regimes where the system was sensitive to intervention, where inhibition of the input signal led to inhibition of the output signal, and where the system was resistant to intervention, that is where the system was robust to input signal inhibition. Moreover, the transition between sensitivity and resistance was governed by a control parameter derived from the relative balance of the activities of three enzymes and drug interventions that target this balance may effect a shift in therapeutic resistance to therapeutic sensitivity.
 
Models such as [20] are powerful approaches to understanding biological mechanisms, since they represent the underlying processes involved. However, they are both data hungry, with hundreds of parameters requiring calibration from (typically) cell line data and/ or estimation based on known system-scale behaviours, and resource intensive to construct. Their development is based on provisional knowledge [21] as not all network associations are known. We propose that the data-driven methods developed here may support process-based model development and investigation. For example, given two sets of time series data from a cell line, where one set represents a control condition and another an intervention condition, our algorithm can construct an interaction network for each of the two data sets. A comparison of the constructed networks would reveal the perturbation manifest in the response of the cell line to the intervention introduced. Any significant differences in the networks would suggest areas of the signalling network most impacted on by that intervention. In turn, this may direct theoretical investigations, such as targetting areas in the signalling network that are highly sensitive to perturbation. More fundamentally, our inference algorithm can support process-based model calibration and validation: time series data from the process-based model should generate similar networks to the cell line data. Finally, it may help interpret model dynamics: by analysing time series data generated from the process-based model in different functioning regimes metalevel associations, not limited to the network topology intrinsic to the process-based model, may be elucidated. For example, an intervention may introduce a marked increase in overall association among a subset of measurables - not necessarily in a localised region of the signalling network topology and this can reveal inform understanding of system-scale responses to that intervention [22].
 
The inference methods presented in this paper are robust and flexible enough to partly incorporate existing knowledge of network structure prior to estimation of model parameters, and this reduces the number of possible networks. It is worth investigating the nature of the factors that contribute to this non-uniqueness in model outcomes. Of note is that for some networks, even when large amounts of data are available, nonuniqueness is still a challenge. In subsequent work we will propose the use of a combinatorial methodological approach to network inference and its application to analysis of pseudo-data generated from real process-based models of biological systems. These combinatorial methods are founded on the work proposed here, and we propose that the methods can serve as fundamental algorithms upon which other variant inference algorithms may be built.
 
References
 






















 
 
 
This article
DOWNLOAD
» XML (95 kB)
» PDF (1,616 kB)
»
Export citation
x
»
Blog this article
   
CONTRIBUTE
» Write a response
» Read other responses
» Publishing with OPG
   
SHARE
» E-mail this article
» Print this article
» Rights and permissions
   
Share
EXPLORE
Related article at
» Pubmed
» DOAJ
» Scholar Google
 
 
 
 
Untitled Document
| More
OMICS Publishing Group is the member of/publishing partner of/source content provider to
All Published content, except where otherwise noted, is licensed under a Creative Commons Attribution License.