Network-based Observability and Controllability Analysis of Dynamical Systems: the NOCAD toolbox

Introduction

In the life sciences, the determination of driver nodes in networks that play a significant role in the emergence or treatment of diseases is an intensively researched field¹. The importance of determining the proper driver nodes, i.e. the ones that ensure physically feasible controllability with the minimum cardinality and energy requirement, in biological networks, or more generally in any dynamical system, is unequivocal, and the amount of research concerning network science has increased rapidly. A detailed study of the control principles in biological networks has already been published². A review about the utilisation of the network science-based determination of driver nodes has also been published that introduced the results of the analysis of the protein-protein interaction (PPI) networks, Caenorhabditis elegans neuronal network, neurochemical rat brain network, Saccharomyces cerevisiae cell cycle networks, Epithelial Mesenchymal Transition (EMT) network, myeloid differentiation regulatory network and Th differentiation network, moreover, the identification of drug targets was also presented³.

The network science-based analysis of dynamical systems has spread rapidly as it provides simple and efficient tools to analyse the structural controllability of any linear or linearised system¹. In terms of controlling the human signalling network, the role of different proteins was also systematically analysed with the toolset of network controllability in 4 to highlight the role of cancer-associated genes. Target control with objective-guided optimisation (TCO) was introduced to control a set of variables (or targets) of interest while the number of drivers and constrained nodes were minimised and maximised, respectively. This method is capable of determining the leading phenotype transitions in biological networks that can be identified as drug targets⁵. In large-scale human liver metabolic networks (HLMN), the driver metabolites have essential functions, moreover, the role of transport reactions and extracellular metabolites in terms of controlling HLMN have revealed the importance of the environment of human liver metabolism with regard to the health of the liver⁶. Using statistical analysis, a subset of critical control nonprotein-coding RNAs (ncRNAs) enriched by human disease can also be determined⁷. In intra-cellular networks, to understand the information flow, a natural control system was utilised and the robustness of such a control was analysed⁸.

The contribution of this paper is to introduce the novel toolbox, NOCAD⁹, and its applicability in the life sciences through the example of the local network of 131 frontal neurons of Caenorhabditis elegans¹⁰. The proposed toolbox is also suitable for the comprehensive analysis of any linear or linearised dynamical systems through their static network representation^11–13. Although in the literature the phrase dynamical network is commonly used, it does not mean that the nodes or connections are temporal but refers to the network of dynamical systems. In the nonlinear case, the methodology needs further clarification because for small nonlinear examples the results can be incorrect¹⁴ and the cardinality of the assigned sensors underestimated¹⁵. As a result, this toolbox deals with only the linear case, nonlinear system-related methods will be implemented later. In the following sections, the representation of linear systems as well as their structural controllability and observability are introduced. Then the theoretical background of the methodology is presented and the implemented functions and measurements introduced through the network of rostral ganglia of C.elegans.

Existing software

Although considerable research has utilised this method¹⁶, a flexible software tool which may be used to support research in this field has yet to be designed. Parallel studies have resulted in a collection of applications, toolboxes, plug-ins and scripts that analyse and determine several structural properties of genes, protein-protein interactions and even social or urban networks. Most of these applications only analyse the structural properties of static networks and just a handful of them utilise these structural properties to draw conclusions concerning the dynamics of the system investigated. As our toolbox belongs to the second group, in the following section, the available applications and programs of this group are elaborated on.

A brief summary of the available tools with expanded functionalities is given in Table 1. Applications or software packages implemented in Python and capable of analysing the controllability and observability of dynamical systems are: graph-control¹⁷ and WDNfinder¹⁸. The advantage of Python-based development lies in its widespread use and the countless methods and packages implemented in this language, including the tools developed for network analysis¹⁹. Although in Python the focus is on developing a broad software package for complex systems analysis, this has yet to be fulfilled and all of the available solutions have limitations. The graph-control toolbox only analyses the impact of network topology on the number of inputs and implements the fast matching algorithm²⁰. Even though WDNfinder only determines the minimum driver node set (MDS) and classifies nodes based on MDS, it is incapable of facilitating extended analysis.

Additionally, the CytoCtrlAnalyser²³ plug-in for Cytoscape²⁵ has been developed, which was implemented in Java and offers graphical user interfaces as well. It evaluates control centrality, control capacity and classifies nodes for biomolecular networks. Furthermore, the Ecological Network Analysis in R software package (enaR) provides some dynamical analysis functions and can generate models to analyse ecological networks in the R environment²⁴. As can be seen, both software packages deal with special kinds of networks. The netctrl program can determine the driver nodes and switchboard dynamics model for any complex network²¹. CONTEST is a MATLAB toolbox which can analyse the dynamics of complex systems, but these dynamics do not cover the structural controllability and observability properties²² of the analysed system. Although the presented software packages ensure the design of a controllable and observable system, they do not provide the opportunity to analyse the designed system exhaustively. These functions are helpful in terms of supporting the work of experts, but are insufficient for the sophisticated analysis of systems.

Methods

In the background of the toolbox the linear systems and their structural controllability and observability properties are stood²⁶. A linear time-invariant (LTI) system is commonly described by its state-space representation that consists of the state equation (Eq. 1) and the output equation (Eq.2).

In the state-space representation, x stands for state variables, u represents the inputs, i.e. the actuators, and y denotes the vector of outputs, i.e. the sensors of the system. Matrices A and B define how state variables and inputs influence changes to the state variables, while matrices C and D define how state variables and inputs influence the outputs, respectively. The cardinality of state variables, inputs and outputs are noted by N, M and K, respectively.

A dynamical system is said to be controllable if it can be driven from any initial state to any desired final state within a finite time with properly selected inputs. Observability is the mathematical dual of controllability. A system is said to be observable if its state can be determined at a given time by a finite set of measured input and output variables. Kalman’s rank criterion²⁷ was used to determine the structural controllability and observability as follows: if the rank of the controllability matrix is equal to the number of state variables, rank(𝓒) = N, then the system is structurally controllable, where 𝓒 = [B, AB, … , A^N–1 B]. Analogously, if the rank of the observability matrix is equal to the number of state variables, rank(𝓞) = N, then the system is structurally observable, where 𝓞 = [C^T, (CA)^T, … , (CA^N–1)T]^T.

To ensure controllability (or observability) using a minimum number of inputs (or outputs), a brute force approach should generate 2^N – 1 configurations of matrix B (or C). To solve this challenging task, the maximum set of disjoint edges is generated by the maximum matching algorithm¹. Two edges are disjointed if they do not share a common starting point or endpoint. The matched nodes are the endpoints of the edges that are a member of the maximum set of disjoint edges, the others are unmatched. Then the unmatched nodes that are generated based on A are the sensor nodes, where outputs should be placed to grant structural observability, while the unmatched nodes generated based on A^T are the driver nodes, where inputs should be placed to grant structural controllability. A^T is also the adjacency matrix of the network representation that is the input of the toolbox. It is very important to note that the result of maximum matching is not unique, and it is possible that the matching is perfect, i.e. no unmatched nodes have resulted. In our implementation, the canonical decomposition of Dulmage-Mendelsohn was utilised to calculate maximum matching²⁸.

For a better understanding, we illustrate the aforementioned definitions by a small example in Figure 1 that contains the command interneurons AVAL, AVAR, AVBL, AVBR, AVDL and AVDR from the frontal neural network of neurons and synapses in C. elegans.

Figure 1. Representation of the command interneurons as a linear dynamical system.

The adjacency matrix of the command interneurons, their network representation and the state equation without assigned input. In this example, due to the symmetric edge pairs between the nodes, the matching is perfect, i.e. all the nodes are matched. In this case, structural controllability and observability can be granted by selecting any node as a driver node and any node as a sensor node.

With the help of the presented Octave- and MATLAB-compatible toolbox, experts can create, analyse and improve any type of dynamical systems. As the structure of the dynamical systems is generally represented by their adjacency matrix and linear dynamical systems can be described by the state-space model that contains the dynamical, input, output and feedthrough matrices, the Octave/MATLAB programming language is a perfect environment to handle these matrices and provide comprehensive functionalities based on them. With the use of NOCAD⁹, experts and researchers can effectively determine the input and output matrices of state-space models, calculate system-specific qualitative measurements (e.g. diameter, relative degree, control centrality and robustness of the system, etc.) and improve the system to satisfy the relative degree-based requirements. The workflow of the toolbox can be seen in Figure 2.

Figure 2. Workflow of the utilisation of the NOCAD toolbox.

The network mapping module provides two methods to create a dynamical system based on the topology of the state variables. The system characterisation module generates more than 49 measures to analyse, classify and characterise the developed system. The improvement and robustness module offers five algorithms to improve the system with additional inputs (observers) as well as outputs (controllers) and can analyse the robustness of the designed system.

Use cases

In this section, the main functionalities of the NOCAD toolbox⁹ are presented through the analysis of the local network of 131 frontal neurons of Caenorhabditis elegans. The first step in the workflow is to create a state-space model based on the adjacency matrix that presents the structural description of the system that, in this case, has the size of 131×131 according to the 131 frontal neurons.

Two methods, path finding and signal sharing are proposed that were implemented to correct the insufficient result of maximum matching. Both methods are modified versions of the maximum matching algorithm. The maximum matching method determined the following 12 neurons to be driver nodes: RMEL, RMER, SIADL, SIADR, SIAVL, SIAVR, SIBDL, SIBDR, SIBVL, SIBVR, SMDDR and URYDR, moreover, determined 12 sensor nodes that correspond to the following neurons: AINL, ASHL, ASIR, ASJR, AWAL, IL2DL, IL2DR, IL2L, SIBDL, URBL, URBR and URYDL. As no critical strongly connected components were present, the results were identical in the case of both the path finding and signal sharing methods.

After utilising the second module of the toolbox, the measures that qualify the whole network with one value are introduced, as presented in Table 2. The network contains 131 neurons and 764 synapses. The density shows that the number of edges is less than a twentieth of the possible maximum, and the diameter of the system, namely the longest shortest path in the network that presents its structure, is 9. The degree variance is 44.3299 which is relatively high given the size of the network, while the Freeman’s centrality is 0.2057. The relative degree of the system is also 4. The Pearson correlation coefficient shows that the in-in, in-out and out-out correlations are slightly assortative in nature, while the out-in correlation is likely to be disassortative. The system is controllable and observable. As no loop is present in the network, the percentage of loops relative to edges is 0%. As 77 symmetrical connections are present between 687 connected node pairs, the percentage of the symmetric edge pairs is 11.2082%.

The second module generates node centrality measures that can reveal structurally important nodes. Since the generated measures can be presented by large tables, they are attached in Excel format to the toolbox⁹. This analysis shows that one of the most important values is the highest degree of the nodes, which belongs to RIAR, an interneuron located in the nerve ring³¹. As Scott’s centrality is a normalised degree, the most important node is once again RIAR. The closeness of node x_i is calculated as the ratio of the number of nodes reachable from x_i to the sum of their distances from x_i . The higher value indicates the more central position of the node, and now RIAL is the most central element. The betweenness centrality shows how many shortest paths intercept the given node. If a node has a high value, then it is a critical node in the structure. The highest value belongs to neurotransmitter RIH that is a serotonin³². The PageRank assigns a percentage value to each node, based on their centrality roles if Markov-chains are modelled. The measure referred to as correlation shows the proportion of the number of edges of neighbours’ and the number of neighbours. This information is useful when determining the assortativity of the system. The control centrality and observe centrality measures determine how many state variables can be influenced or observed by the nodes.

The determined driver and sensor nodes can be classified into four groups³³. According to these groups, four phenomena can provide driver or sensor nodes. Firstly, source nodes when the node has no incoming edges, thus, a dedicated input is needed. Secondly, dilation, when the generated set of child nodes has higher cardinality than the number of parent nodes. A distinction is made between internal dilation and external dilation, in the former the child node is not a leaf, i.e. it has children, while in the latter the child is a leaf node, i.e. it has no children. The last type is the inaccessible nodes when the node has an incoming edge and no dilation is present, but the node is not reachable by a directed path from any of the inputs. These types are important properties, e.g. the existence of dilation or inaccessibility is detrimental to complete structural controllability³. The controlling and observing matrices are sparse matrices as only the columns of drivers and sensors contain nonzero values. The values show the number of derivations necessary to influence or observe a state variable in the system. Next, the similarity of the driver and sensor nodes is presented. This similarity is based on how similar the set of nodes is, which can be reached for driving or observing. Furthermore, the necessary derivation to influence or observe them is also part of the comparison. R_c and R_o are the simple reachability matrices. They show which nodes can be controlled or observed by a given node in its structural meaning, i.e. the existence of a directed path between the nodes is shown. In R_c, the i^th column shows which nodes can control node i. From the other viewpoint, elements in row i highlight those nodes which can be controlled by node i. It is very important that R_c is only a reachability matrix, the structural controllability of the reachable nodes is not granted by a node that can reach them, but in some cases the structural controllability problem can be reduced to a reachability problem³⁴. The R_o matrix can be interpreted analogously with regard to observability.

Finally, measures of edge centrality are generated by the system characterisation module. The betweenness has the same meaning as in the case of nodes, that is, it yields the number of shortest paths that intercept the edge³⁵. From this perspective, the most critical synapsis is the one between the command interneuron AVAL and amphid ADLL with a value of 640.5833. The endpoint similarity shows how similar the influenced and observed sets of the state variables with regard to the endpoints of edges are. This metric has a high value if the edge is part of a cycle or creates a bridge in the network. As no bridges are present in this network, only cycles can be recognised by this measure. The edge similarity shows how similar the roles of edges are, and it allows redundancies, to be located.

For the demonstration of the last module, four plus one methods were applied to the neural network of C. elegans. The set covering-based grassroot method (SetCovGr) optimises the placement of driver nodes and sensor nodes to provide an initially demanded relative degree, but this method does not take into account the original input and output configurations also, thus, structural controllability and observability is not granted also. The other four methods grant controllability and observability by expanding the minimal configurations. They are the centrality measures-based (CentMeas) retrofit, set covering-based retrofit (SetCovRet), modified Clustering Large Applications based on Simulated Annealing (mCLASA) and Geodesic Distance-based Fuzzy c-Medoid Clustering with Simulated Annealing algorithm (GDFCMSA) methods¹³. These methods were utilised with the following parameters: the required relative degree was set at 2, while the alpha parameter of the cost function was set at 0.5¹³. The results can be seen in Table 3. The number of assigned driver nodes varies significantly when different methods are applied. The centrality measures-based method assigned the most driver nodes to the system. Thus, this method results in the smallest cost, but the difference is irrelevant, most of the methods resulted in a cost of 1.5. The increase of the number of the driver nodes decreases the mean relative degree, which is the lowest in the case of the centrality measures-based method.

The robustness of the configuration was also analysed. In each scenario, a node was removed from the network. Using the leave-one-out strategy, the network with the altered configuration remains controllable in 115 scenarios. As for the sensor nodes, the difference is not as significant between the methods as in the case of the driver nodes. Critical nodes were also generated. A node is critical if the system becomes uncontrollable or unobservable if the node is removed. The determined critical nodes and the names of selected driver and sensor nodes can be found in the Excel file attached to the toolbox.

Conclusions

Although numerous papers have utilised the network-based determination of driver and sensor nodes, a flexible toolbox that may be used to support the analysis has yet to be designed. To fill this gap, in this article the Octave- and MATLAB-compatible NOCAD toolbox⁹ was proposed to support the network-based controllability and observability analysis of dynamical systems, and through the analysis of the neural network of C.elegans, the applicability of the toolbox in the life sciences was presented. The toolbox offers two methods to design a structurally controllable and observable system based on the adjacency matrix (A^T). The designed system can be analysed by 49 qualitative measures both from structural and dynamical points of view. The toolbox serves five methods to improve the designed system by adding new inputs and outputs to it, thus, its relative degree can be decreased. Then the robustness of the individual designs can also be evaluated. The modular structure of the toolbox supports the facile improvement of the modules by adding new functions and the toolbox can be extended by new modules as well. Even though the modules are built on each other, most of their functions can also be used independently from each other.

Although our goal in this paper is to draw the attention of researchers of life sciences to the services provided by the NOCAD toolbox, it can be utilised in practice in various fields of sciences as well, for example, it enables social networks to be controlled in the economy, transaction networks to be analysed in finance or dynamical systems to be designed in engineering.

Data availability

All data underlying the results are available as part of the article and no additional source data are required.

Software availability

Source code available from: https://github.com/abonyilab/NOCAD.

Archived source code at time of publication: https://doi.org/10.5281/zenodo.2656674⁹

License: GNU General Public License v3.0