Analysing Gene Networks with PDP Systems. Arabidopsis thaliana, a Case Study

Analysing Gene Networks with PDP Systems. Arabidopsis thaliana, a Case Study Luis Valencia–Cabrera1, Manuel Garc´ıa–Quismondo1, Mario J. P´erez–Jim´enez1, Yansen Su2, Hui Yu2, and Linqiang Pan2 1Research Group on Natural Computing Dpt. of Computer Science and Artificial Intelligence, University of Sevilla Avda. Reina Mercedes s/n. 41012 Sevilla, Spain {lvalencia,mgarciaquismondo,marper}@us.es 2Department of Control Science and Engineering, Huazhong University of Science and Technology 1037 Luoyu Road, Wuhan, China [email protected], [email protected], [email protected] Summary. Gene Regulatory Networks (GRNs) are a useful tool for biologists to understand the interactions among genes in living organisms. A special kind of GRNs known as Logic Networks (LNs) has been recently introduced. These networks consider that the state of one or more genes can influence another one. In a previous work, we proposed a Membrane Computing model which simulates the dynamics of LNs by drawing on the improved LAPP algorithm. In this paper we provide a case study for our LN model on a network which regulates the circadian rhythms of long–term studied plant Arabidopsis thaliana. We outline the software tools employed and propose a methodology for analysing LNs on our Membrane Computing model. At the end of the paper, some conclusions and future work are included. Keywords: Bioinformatics, Genetics, Gene networks, Membrane Computing, MeCoSim, Software engineering, Modelling, LAPP, Logic networks 1 Introduction Since its very beginning, Membrane Computing [13] has been employed as a modelling framework for biochemical phenomena. Although the current landscape is more focused on metabolite–oriented dynamics, gene regulatory networks (GRNs) have also been modelled by means of P systems as part of this framework. In a previous work, we followed this line of research by proposing a Membrane Computing model for a specific type of gene networks known as Logic Networks (LN ) [16]. This model describes a P systems family known as LN Dynamic P systems (LN 258 L. Valencia–Cabrera et al. DP systems), within the framework of PDP systems [11]. LNs are a specific type of GRNs in which the combination of states of several genes, rather than the single state of any of them, influence another one. Bowers et al. [2] proposed a methodology for the construction of logic networks out of statistical data, known as Logical Analysis of Phylogenesis Profiles (LAPP). In our model, these combinations are limited to at most two genes affecting a third one. The model, in conjunction with DCBA algorithm[4], intends to capture the behaviour of the Improved LAPP Method introduced by Wang et al. [18]. In their work, they propose a case study on a gene network associated to Arabidopsis thaliana’s flowering process. We intend to reproduce this case study by using our Membrane Computing model. We also include a guide for generating custom simulators on MeCoSim for LN DP systems, depicting a step–by–step guide on MeCoSim tool [12]. Finally, the data employed in this case study is provided as an appendix, thus easing cross–checking of results. This paper is structured as follows. Section 1.1 introduces the Logic Network to be studied, a GRN associated to the flowering process of Arabidposis thaliana. Section 2 outlines the LN DP system model presented in [16], in order to make the current work self contained, as it is used to analyse our case study. Section 3 consists of a guide to simulate LNs from scratch on MeCoSim [14, 12]. This guide complements the simulation methodology described in [16]. Section 4 describes a case study on a real–world logic network on Arabidopsis Thaliana, in order to experimentally verify the behaviour of the model on complex gene networks. Finally, section 5 lists the conclusions obtained and proposes some open problems. 1.1 A Logic Network on Arabidopsis thaliana flowering processes Arabidopsis is a long–day plant. Zhang and Zuo [19] stated that long–day conditions can promote reproductive growth and induce early flowering. However, short–day conditions can promote vegetative growth and induce late flowering or even no flowering. To understand the intrinsic mechanisms of Arabidopsis flowering in different lighting conditions, it is required to compare the relationships of related genes. In the latest ten years, much work has been reported in the field about A. thaliana flowering. Imaizumi et al. [8] found that FKF1 is a blue light receptor which regulates flowering. Later, they also showed that FKF1 together with Flavin–Binding and Kelch Repeat degrade Cycling Dof Factor1 (CDF1) to eventually control carbon monoxyde [7]. In the same year, Abe et al. [1] found that Flowering Locus T (FT) together with FD activate Apetala1 (AP1) to initiate floral development and promote floral transition at the shoot apex. Previous work deal only with one or few genes related to flowering. However, the networks considered in this work focus on the relationships among a large number of genes systematically. Bowers et al. [2] proposed the Logic Analysis of Phylogenetic Profiles (LAPP) [2]. This method helps researchers to know biological functions of some genes or proteins on the basis of phylogenetic profiles, which has been developed both on theory and application ([3, 20, 17]). For example, Wang et al. [17] Analysing Gene Networks with PDP Systems 259 developed the improved LAPP method, and reversely constructed a logic network of sixteen genes in shoot for Arabidopsis under salt stimuli. 2 Description of the model This section summarizes both the P system family and the model (i.e., initial configuration and rule patterns) employed in this case study. For a detailed description of the model, see [16]. 2.1 A family of P systems based on Logic Networks The model depicted here is a P system of a family known as Logic Network Dynamic P systems (LN DP systems). An LN DP system is described within an expansion of Population Dynamics P systems (PDP systems) [16]. An LN DP system ΠLN of degree (q, m) with q, m ≥1, taking T≥1 time units, is a tuple ΠLN = (G, Γ, Σ, T, RE, µ, R, {fr,j :r∈R, 1≤j≤m}, {Mij : 1 ≤i≤q, 1≤j≤m},{Mj: 1 ≤j≤m}) where: •(G,Γ,Σ,T,RE,µ,R,{fr,j :r∈R, 1≤j≤m},{Mij : 1 ≤i≤q, 1≤j≤m}) is a PDP system. • {fr,j =1: r∈R, 1≤j≤m}. •For each j(1 ≤j≤m), Mjare multisets over Γ, describing the objects initially placed in the menvironments ej. In this paper, in the description of an LN PDP System, functions fr,j are omitted. They are all equal to 1, so it is not necessary to make them explicit. 2.2 The model Here the model for the family of Logic Network Dynamic P systems is outlined. This model covers any possible P system in this family, so the multisets, rules, etc. depend on the P system which represent each specific instance of a logic network. The definition of the general model requires the use of parameters in our constructs, as explained at the end of this subsection. Let LN be a logic network. Let ng,nu,nb be the number of genes, unary and binary interactions of LN, respectively. Let n=ng +nu +nb. The model consists of the following PDP system of degree (1, n), ΠLN = (G, Γ, Σ, T, RE, µ, R, {Mij : 0 ≤i≤q−1,1≤j≤m},{Mj: 1 ≤j≤m}) where: 260 L. Valencia–Cabrera et al. •Gis a directed graph containing a node (environment) for each gene, unary or binary interaction, following this order. •In the alphabet Γ, we represent gene states, interaction types, contribution weights and targets. Γ={ai, bi, ci: 0 ≤i≤1} ∪ {go, d0} ∪ {unopj, binopj: 1 ≤j≤4} ∪ {auxDesti,gj,1,k : 0 ≤i≤1,1≤j≤ng, 1≤k≤nb +nu} ∪ {desti,gj,1,tk,1+ng : 0 ≤i≤1,1≤j≤ng, 1≤k≤nb} ∪ {desti,gj,1,untk−nb,1+ng+nb : 0 ≤i≤1,1≤j≤ng, nb + 1 ≤k≤nb +nu} ∪ {etk,4∗i+(1−i)∗(1−tk,4),tk,1+ng : 0 ≤i≤1,1≤k≤nb} ∪ {etk,6∗i+(1−i)∗(1−tk,6),tk,1+ng : 0 ≤i≤1,1≤k≤nb} ∪ {euntk−nb,4∗i+(1−i)∗(1−untk−nb,4),untk−nb,1+ng+nb : 0≤i≤1, nb + 1 ≤k≤nb +nu} ∪ {eFtk,8∗i+(1−i)∗(1−tk,8),tk,1+ng : 0 ≤i≤1,1≤k≤nb} ∪ {eFi,(untk,1+ng+nb): 0 ≤i≤1,1≤k≤nu} ∪ {clockj: 0 ≤j≤cc + 3} •The environment alphabet is Σ=Γ\ {d0} •Each cycle to evolve from a real network configuration to the next one involves 15 computational steps, so T= 15 ·Cycles, where Cycles is the number of cycles to simulate. •µ= [ ]1is the membrane structure. •The initial multisets are: –Mgk,1={a1gk,3, a01−gk,3,go : 1 ≤k≤ng}. –Mng+ti,1={binopti,2: 1 ≤i≤nb}. –Mng+nb+unti,1={unopunti,2: 1 ≤i≤nu}. •The rules of Rand REto apply are showed below. They are put together to follow the sequential order of execution. Environment rules start with re and skeleton rules start with rs. –rs1,i ≡go ai[ ]1−−→ cibimax∗ib0thresholdclock0[ ]1: 0 ≤i≤1 – For each source gene environment: re2,i,j,k ≡(ci−−→{auxDesti,gj,1,k :{1≤k≤nb +nu}})gj,1 : 0 ≤i≤1,1≤j≤ng re3,i,j,k ≡(auxDesti,gj,1,k −−→ desti,gj,1,tk,1+ng )gj,1 : 0 ≤i≤1,1≤j≤ng, 1≤k≤nb re4,i,j,k ≡(auxDesti,gj,1,k −−→ desti,gj,1,untk−nb,1+ng+nb)gj,1 : 0 ≤i≤1,1≤j≤ng, nb + 1 ≤k≤nb +nu re5,i,k ≡(desti,tk,3,tk,1+ng −−→ etk,4∗i+(1−i)∗(1−tk,4),tk,1+ng )tk,3 : 0 ≤i≤1,1≤k≤nb re6,i,k ≡(desti,tk,5,tk,1+ng −−→ etk,6∗i+(1−i)∗(1−tk,6),tk,1+ng )tk,5 : 0 ≤i≤1,1≤k≤nb re7,i,k ≡(desti,untk−nb,3,untk−nb,1+ng+nb −−→ euntk−nb,4∗i+(1−i)∗(1−untk−nb,4),untk−nb,1+ng+nb)untk−nb,3 : 0 ≤i≤1, nb + 1 ≤k≤nb +nu Analysing Gene Networks with PDP Systems 261 re8,i,k ≡( )tk,1+ng(ei,tk,1+ng)tk,3−−→(ai)tk,1+ng ( )tk,3 : 0 ≤i≤1,1≤k≤nb re9,i,k ≡( )tk,1+ng(ei,tk,1+ng)tk,5−−→(ai)tk,1+ng ( )tk,5 : 0 ≤i≤1,1≤k≤nb re10,i,k ≡( )untk−nb,1+ng+nb(ei,untk−nb,1+ng+nb)untk−nb,3−−→ (ai)untk−nb,1+ng+nb( )untk−nb,3 : 0 ≤i≤1, nb + 1 ≤k≤nb +nu – Evaluation of the result of the interactions (1/2). rs11 ≡binop1a02[ ]1−−→ binop1c0[ ]1 rs12 ≡binop1a12[ ]1−−→ binop1c1[ ]1 rs13 ≡binop1a1a0[ ]1−−→ binop1c1[ ]1 rs14 ≡binop2a12[ ]1−−→ binop2c1[ ]1 rs15 ≡binop2a02[ ]1−−→ binop2c0[ ]1 rs16 ≡binop2a1a0[ ]1−−→ binop2c0[ ]1 rs17 ≡binop3a12[ ]1−−→ binop3c0[ ]1 rs18 ≡binop3a02[ ]1−−→ binop3c0[ ]1 rs19 ≡binop3a1a0[ ]1−−→ binop3c1[ ]1 rs20,i ≡unop1ai[ ]1−−→ unop1ci[ ]1: 0 ≤i≤1 rs21,i ≡unop2ai[ ]1−−→ unop2ci−1[ ]1: 0 ≤i≤1 rs22,i ≡unop3ai[ ]1−−→ unop3cii[ ]1: 0 ≤i≤1 rs23,i ≡unop4ai[ ]1−−→ unop4c1−ii[ ]1: 0 ≤i≤1 – Evaluation of the result of the interactions (2/2). re24,i,k ≡(ci)tk,1+ng( )tk,7−−→ ( )tk,1+ng(eFtk,8∗i+(1−i)∗(1−tk,8),tk,1+ng )tk,7 : 0 ≤i≤1,1≤k≤nb re25,i,k ≡(ci)untk,1+ng+nb( )untk,5−−→ ( )untk,1+ng+nb(eFi,(untk,1+ng+nb))untk,5 : 0 ≤i≤1,1≤k≤nu – Calculation of contributions. rs26,i,k ≡eFi,(tk,1+ng)[ ]1−−→ bitk,9[ ]1: 0 ≤i≤1,1≤k≤nb rs27,i,k ≡eFi,(untk,1+ng+nb)[ ]1−−→ biuntk,6[ ]1: 0 ≤i≤1,1≤k≤nu – Elimination of different–signed contributions. rs28 ≡b1b0[ ]1−−→[ ]1 rs29,i ≡clocki−1[ ]1−−→ clocki[ ]1: 1 ≤i≤cc + 3 – Calculation of the next gene state. rs30 ≡b0[ ]− 1−−→[d0]− 1 rs31 ≡b1[ ]− 1−−→[]− 1 rs32,i,j,k ≡desti,j,tk,1+ng[ ]− 1−−→[ ]− 1: 0 ≤i≤1,1≤j≤ng, 1≤k≤nb rs33,i,j,k ≡desti,j,untk−nb,1+ng+nb[ ]− 1−−→[ ]− 1 : 0 ≤i≤1,1≤j≤ng, nb + 1 ≤k≤nb +nu rs34 ≡[d0]− 1−−→[ ]+ 1 rs35 ≡clockcc+3[ ]+ 1−−→ go a0[ ]0 1 rs36 ≡clockcc+3[ ]− 1−−→ go a1[ ]0 1 In this section, only input parameters are described. This way, details about the model dynamics are left aside. These parameters are described in table 1. 262 L. Valencia–Cabrera et al. Table 1. Parameters Parameter Description General parameters for the system ng Number of genes in the network nb Number of binary interactions nu Number of unary interactions threshold Maximum strength for an interaction cc Clock control Gene configuration parameters gi,1Gene number (id) gi,3Initial state of the gene Binary interactions parameters ti,1Binary interaction number (id) ti,2Interaction type (or: 1, and: 2, xor: 3) ti,31st source gene number (id) ti,41st source gene contribution (positive: 1, negative: 0) ti,52nd source gene number (id) ti,62nd source gene contribution (positive: 1, negative: 0) ti,7Destination gene number (id) ti,8Influence over destination gene (positive: 1, negative: 0) ti,9Strength of the destination Unary interactions parameters unti,1Unary interaction number (id) unti,2Interaction type (strong promotion: 1, inhibition: 2; weak ones: 3, 4) unti,3Source gene number (id) unti,4Source gene contribution (positive, negative) unti,5Destination gene number (id) unti,6Influence over destination gene (positive, negative) 2.3 Model output The state of the network is encoded as the multiplicity of objects a1and a0in each gene environment. The presence of objects a1inside a gene environment represents that its gene is active (a0for inactive). Due to the nature of the system, membrane genes cannot have objects a1and a0simultaneously. Therefore, to know the final state of the network, it suffices to identify which environments contain object a1 and which ones a0at configuration T. Analysing Gene Networks with PDP Systems 263 3 Modelling and simulation on MeCoSim This section explains some relevant issues concerning the software environment, putting the focus on the needed changes in P-Lingua framework and the configuration of MeCoSim. The P–Lingua definition used to analyse the PDP model adheres to P–Lingua version 4 standard, available at [10]. 3.1 Custom interface in MeCoSim In our previous work [16], a P–Lingua model for the family of logic networks based on PDP systems has been extensively described. This model contains a number of parameters representing relevant information about each specific scenario. Thus, although a general model has been presented, a mechanism to ease the task of introducing the specific data for each scenario is needed. This task is performed through the software environment provided by MeCoSim [14, 12]. MeCoSim permits the definition of a custom visual simulator. This simulator includes an interface with the needed inputs, outputs, and a way to translate the input data into parameters for the model. The simulation engine is provided by pLinguaCore, available at [10]. The most relevant facts of this process are listed below. Definition of a custom visual simulator for Logic Networks Here, the process for defining a custom simulator based on MeCoSim is provided. This process is very simple, and consists of the following steps: Configuration file: The first step is to define a spreadsheet file containing the configuration for the definition of visual tabs, input tables, output tables and charts, and the mechanism to generate both model parameters from input tables and outputs from the simulation results. The contents of the simulation parameters tab in the file is shown in figure 1. The configuration file is available by contacting the authors. Fig. 1. MeCoSim configuration file. Simulation params Loading configuration file on MeCoSim: That file is loaded through the main window in MeCoSim by clicking the “Load config file” button, choosing the file, selecting “Update all information” option and pressing “Update config info” button. After these steps, the configuration file is loaded, so the custom simulator is ready to use. Finally, the message “The Application has been successfully initialized” is prompted in MeCoSim main display. 264 L. Valencia–Cabrera et al. Running custom simulator: The newly configured simulator is ready to use by selecting “Gene network” application and pressing “Run Application” button. Then, the custom interface is visualized, enabling the user to load the model (.pli file) and enter the input data for a specific scenario, as shown in figure 2. Fig. 2. MeCoSim window 3.2 Simulation methodology In [16], we describe a methodology to simulate LN DP systems on MeCoSim. This methodology can be summarized in the following steps: •Load the model specification by clicking on Model > Set model. •Fill in the input tables in tab Input. Optionally, it is possible to save this data by clicking on Scenario > Save. This data can be loaded later by clicking on Scenario > Open. •Set the number of steps on Simulation > Number of steps. •Click on Simulation > Simulate!. •Visualize the results in tab Output. A toy example on a 3–gene logic network is provided in [16]. This network is taken from [15]. In this network, interactions have no associated weights. Hence, we presume all interactions to have the same weight (say 100). Although interaction scoring based on Pearson correlation coefficient is a rather widespread metric for measuring gene interaction strength [9], there is little literature on LNs, thus making it hard to find LN toy examples. Analysing Gene Networks with PDP Systems 265 4 A case study on Arabidopsis thaliana In order to experimentally verify our model, we have tested our algorithm by using a logic network which regulates flowering processes associated to Arabidopsis thaliana on a long day scenario. This relatively large network integrates gene interaction samples from NCBI/EBI database [5]. This logic network has been constructed according to the procedure described by Bowers et al. [2]. A. thaliana is a species widely used in genetic and protein interaction networks. The total number of genes in the network is 29, whereas the total number of interactions is 99. These interactions consist of 23 unary interactions and 76 binary interactions. We notice that only a few different types of all possible interactions are present in this network. In the case of unary interactions only strong promotions and strong inhibitions are present. When it comes to binary interactions, only AND– like and OR–like interactions are present. As regards to the distribution of the present interactions, the vast majority of them are AND–like interactions with both inputs in non–negated form (that is, G0 j=Gjand G0 k=Gk), as well as non–negated result (G0 l=Gl). Gene network data is provided as an appendix in section 7. Specifically, gene initial states are reflected in table 5. Unary gene interactions are reflected in figure 6. Similarly, binary gene interactions are reflected in figures 7, and 8. Figure 3 displays the MeCoSim input tables used in this case study. Eventually, we have simulated the corresponding P system for the A. thaliana network entered. The improved LAPP method (as presented in Wang et al. [18]) has been run for 30 steps on this data. Similarly, the LN DP model has been simulated for 30 cycles. As each cycle in an LN DP system consists of 15 computation steps, the total number of steps simulated in the model is 30 ×15 = 450. The results (see figure 4) match the ones obtained from the execution of the improved LAPP method on the same input data. Therefore, it is verified that, on this gene network and scenario, the P system model behaviour is analogous to that from the improved LAPP method.