Rule-based multi-scale simulation for drug effect pathway analysis

Computational model for drug response simulation

The simulation model of biological networks could be differed according to the specification of nodes in a system. Nodes that are specified by discrete values are simulated by qualitative method such as Boolean network or Petri-net model. Others nodes that have continuous values are simulated by quantitative method such as ordinary differential equation (ODE) model.

Jack et al. conducted qualitative investigation of cellular responses using Boolean network [10]. They investigated the response of biological systems by the effect of chemicals at molecular level and tissue level. Jin et al. developed enhanced Petri-net model to investigate a genetic mechanism for effects of drug combinations [11]. They integrated data on large-scale biological network and made Petri-net model from the biological network. Although their model is large enough to investigate drug actions, it predicts only cellular response of pairwise drug combinations by Petri-net simulation result.

To obtain accurate results in analysing drug action by biological networks, constructing an accurate biological network is needed. Nelander et al. made a model that derives network models from molecular profiles of perturbed cellular systems [12]. They predicted results of multiple genetic alterations by ODE model.

Previous research focused on constructing models that investigate the biological response of the drug effect. However, previous research did not investigate at whole-body level.

Multi-scale modelling

Type 2 diabetes (T2D) is caused by a problem in the way body makes or uses insulin [13]. Usually when we get T2D, our fat, liver, and muscle cells do not respond properly to insulin. T2D cannot be understood by investigating only one organ, which produces insulin for drug discovery processing. It occurs when several organs' interaction does not work properly. So, we need to investigate systematically organs' interactions. There are many complex diseases that have multi causalities, such as schizophrenia, stroke, etc.

To investigate these complex diseases, we need to model biological systems as multi-level system that contains molecular level, cell/tissue-level and organ-level. Cell-level are affect by biological processes at the molecular level and biological processes at the cell/tissue level influence dynamics at organ level as Figure 1. This hierarchical organization and the causalities between different levels are characteristics of biological systems [14, 15]. As we mentioned, we are lack of understanding about tissue level and organ level. Hence, it is difficult to make a whole-body model mathematically, like in ODE. Qualitative modelling could be used to overcome this difficulty. Among qualitative modelling methods, we will utilize the context of rule-based modelling to model a whole-body level.

Simpler form of rule-based modelling [16, 17] makes it possible to converts rules based on various sources of knowledge, such as literature and computational models of different formalisms. And this formalism can describe more diverse states of a component in biological system than the Boolean network model.

We previously proposed rule-based multi-level simulation model [18]. In previous work, we showed the possibility of our model to simulate the multi-compound drug effect on human body, which is multi-level system. In this work, we improved previous model to have better performance.

Drug effect pathway

Our simulation platform shows drug effect pathways, where the drug actions propagates from the target to the final nodes for the therapy that can change disease state to normal state. Figure 2 demonstrates the metformin effect pathway on T2D model. It shows integrated drug effect pathway of metformin in the rules [20].We simulated effect of single drug at each simulation step among twenty-two drugs. According to simulation results, all twenty-two drugs repress the state of glucose in the circulation system. It could also investigate effects of multi-drugs through a simulation.

According to the metformin effect pathway (Figure 2), metformin targeted PRKAB1 enzyme in hepatocyte. The metformin activate insulin secretion through GLUT4 and inhibits gluconeogenesis, which generate glucose, by inhibiting PEPCK and Glucose 6-phosphatase. This result shows corresponding mechanisms of metformin, which previously studied [23]. Drugs (i.e. pioglitazone, rosiglitazone, troglitazone) targeting muscle cells inhibit insulin resistance and glucose level. Drugs targeting pancreas cells (i.e. exenatide, glyburide etc.) help insulin secretion and regulation of glucose metabolism.

Combination drug

In a case of metformin and rosiglitazone (Figure 3), metformin makes the state of gene PRKAB1 in hepatocyte cell up and rosiglitazone makes the state of PPARγ in adipocyte up. The result of metformin action reduces glucose level in circulation and the result of rosiglitazone action reduces glucose level by different pathways. In this way, we can assist to investigate the mechanism of combination drugs by our simulation result.

Components and rules

The description of a biological system in this work consists of a collection of components and rules. A component in a biological system represents an organ, a cell, or a drug and so on (Figure 4). A rule is composed of left and right side to describe the action of rule:

$\left(\mathsf{\text{Component}}\mathsf{\text{type}},\mathsf{\text{Component}}\mathsf{\text{name}},\mathsf{\text{Attribute}}\mathsf{\text{type}},\mathsf{\text{Attribute}}\mathsf{\text{name}},\mathsf{\text{Condition}}\right)\to \left(\mathsf{\text{Component}}\mathsf{\text{type}},\mathsf{\text{Component}}\mathsf{\text{name}},\mathsf{\text{Attribute}}\mathsf{\text{type}},\mathsf{\text{Attribute}}\mathsf{\text{name}},\mathsf{\text{Action}}\right)$

The condition of a rule signifies the change of state value of an attribute.

The description of condition 'up' is:

The description of condition 'down' is:

The description of action 'up' is:

The description of action 'down' is:

Gene AGRP, which is in hepatocyte cell, change state from 'normal', 'high' to 'high', or from 'low' to 'high' when gene PPKAB1, which is in hepatocyte cell, has changed state from 'normal', 'high' to 'high', or from 'low' to 'normal'.

Rule extraction strategy

One of the main issues of system development is collecting rules of the system in a whole-body scale for the simulation. For this purpose, we have extracted rules from disparate information sources to increase coverage: pathway database, ordinary differential equation (ODE) model and Petri-net.

A rule curated from experts was described based on dynamic status rather than the static. It is based on the assumption that the body system changes its status through perturbation. If a certain component of cellular system does not alter, other components connected to it cannot be affected permanently. Thus, curators have extracted the rules from the context of the data that describes perturbed status of a certain component. A well-curated signalling pathway database, pathway interaction database (PID) provides valuable information of diverse signalling protein interactions. Regulations of proteins in each pathway were annotated explicitly, so it could be converted to our rule formulae straightforwardly. For the analysis of Type 2 diabetes (T2D), related pathways like insulin pathway and insulin-mediated glucose transport pathway in PID were considered to be extracted as a rule for the simulation. This database also provides integrated resources from external databases, such as BIOCARTA, REACTOME. Therefore, we could also collect other insulin pathway-related rules of external databases through these databases.

Rule conversion from different modelling formalisms

An ODE model is composed of several differential equations describing relationship between variables in the model. Some of the variables have clear relationship (i.e., activation or inhibition) in, for instance, first-order linear differential equations. However, some of equations of nonlinear form make the relationship of variables difficult to determine.

We have examined whether the effect of a certain variable to others shows a monotonic changes. Taking an advantage of rule-based simulation, we have only focused on the direction rather than the quantity. We were able to extract rules of variables if they affect each other monotonically. We also identified the nonmonotonical relationship of variables if the range of changing pattern of effect could be determined.

Conventional Petri-net has the graphical structure containing a transition and a place as an element of the graph. A transition element means the changing action, not the actual component, that mediates the interconversion of place elements of the Petri-net through its action, and a place element is a component of the body system that takes the change of its status. Therefore, based on the context of transition elements, we were able to extract rules from places elements (i.e., place elements will be located in the left- and right-hand-side of a rule with its perturbed status determined by the context of transition elements.)

The overview of simulation

This system has three kinds of inputs, such as a component, a rule, and a disease file. The component file contains information about components such as component types, component names, attribute types and attribute names. The system constructs components according to the component file (Figure 5). The disease file contains disease states of components. The system initiates the state of components by a disease file to make the system as a patient model. After all components in the system are initiated, the system simulate according to the rules in the rule file. The result file will be produced after simulation is over. The result file has all state changes of every component during the simulation. Each component has the threshold for changing its state. Each threshold of the components is different depending on the component types.

Simulation algorithm

This system uses an asynchronous updating simulation method to observe drug responses and drug effect pathways to implement stochastic rule firing. Before running simulation, the system is initialized by T2D patient model and T2D related drugs. Using T2D patient model, the states of components, which are related to T2D, in the system are changed to abnormal. Then, drugs are injected to the simulation system. The action of drug injection means the start of the simulation.

Algorithm.1 and 2 show the simulation running process. R means all rules in the system. C means all components in the system. AR(t) is a set that has all executable rules in the specific time t. AC(t) is a set which has all active components at time t. And RFS(c) is the rule firing score of a component c. After a rule execution, the RFS of a component, which is affected by the rule effect score, which is RS(r), will be changed. TH(c) is the rule execution threshold of a component. If RFS(c) is greater than TH(c), the state of component c is updated by the rule.

To model differences in the speed of signal propagation in the body, we used a random-order asynchronous updating algorithm. Executable rules are selected through Algorithm 1. Next, one rule is randomly selected from executable rules and executed. Then, the states of components are updated if condition of the randomly selected rule is satisfied. Then, all rules activated by the randomly selected rule are added to AR for the next iteration step. We repeated the above simulation 100 times as a Monte Carlo approach for finding drug effect pathways.

R = {x|x is a rule}

C = {x|x is a component}

AR(t) = {x|x is an active rule at time t}

AC(t) = {x|x is a component at time t}

RFS(c) = (rule firing score of a component)

TH(c) = (rule execution threshold of a component)

RS(r) = (rule effect score)

Algorithm.1: RunSimulation( )

Algorithm.2: UpdateComponentSet(AC(t),R _i )

Algorithm.2 shows how to update the states of the components in AC when $R_i$ is executed at time t. First, save all components affected by $R_i$ at $A{C}_i$. Then, repeat update of RFS of $C_j$ and check RFS is greater than threshold of the component $C_j$, which is TH, until AC has not any element. If the RFS is greater than the threshold, the state of the component is updated.

Algorithm.2 shows how to update the states of the components in AC when $R_i$ is executed at time t. First, save all components affected by $R_i$ at $A{C}_i$. Then, repeat update of RFS of $C_j$ and check RFS is greater than threshold of the component $C_j$, which is TH, until AC has not any element. If the RFS is greater than the threshold, the state of the component is updated.