Simulation with ODE
Using the topology of the basic apoptosis network (Figure 1), constructed from the literature, we derived a computational model by generating rate equations for the network. The complete set of differential equations for the basic network N1 is given in Table 2. Concentration of each molecule is represented by enclosing the molecule name with square brackets (e.g. [Casp8]). Reaction rate constants are represented as rc1, rc2 etc. For example, based on the basic network N1 in Figure 1, the differential equations for the change of concentration of AKT and p53 were derived as:
Table 2  Differential equations for the basic apoptosis network N1 (Figure 1)
S.N.  Rate of Change  Differential equations
1  d  d t   Casp  6  - [Casp6]·rc11 + [Casp9]·rc12
2  d  d t   CDK  1  +[Wee1]·rc13 - [CDK1]·rc14
3  d  d t   Wee  1  - [Wee1]·rc13
4  d  d t   JNK  - [JNK]·rc15 + [EGFR]·rc17
5  d  d t   BID  - [BID]·rc16 - [BID]·rc25
6  d  d t   EGFR  - [EGFR]·rc17 - [EGFR]·rc32;
7  d  d t   PUMA  - [PUMA]·rc19 + [p53]·rc27
8  d  d t   Casp  3  - [Casp3]·rc1 + [XIAP]·rc2 + [SMAC]·rc4 + [Casp8]·rc7 - [Casp3]·rc8 +[Casp9]·rc9 +[Casp6]·rc11 + [MOMP]·rc30
9  d  d t   Casp9  - [Casp9]·rc9 - [Casp9]·rc12 + [PUMA]·rc19 - [Casp9]·rc24 + [BCL2]·rc26
10  d  d t   HER2  - [HER2]·rc20 - [HER2]·rc31
11  d  d t   Casp8  - [Casp8]·rc3 + [RIP1]·rc6 - [Casp8]·rc7 + [Casp3]·rc8 + [BID]·rc25 +[EGFR]·rc
12  d  d t   Apop  +[Casp3]·rc1
13  d  d t   XIAP  - [XIAP]·rc2 + [SMAC]·rc21 + [Casp9]·rc24 + [HER2]·rc31
14  d  d t   SMAC  +[Casp8]·rc3 - [SMAC]·rc4 + [BCL2]·rc5 - [SMAC]·rc21 + [MOMP]·rc29
15  d  d t   BCL2  - [BCL2]·rc5 + [p53]·rc10 + [CDK1]·rc14 + [JNK]·rc15 +[BID]·rc16 +[STAT3]·rc22 - [BCL2]·rc26
16  d  d t   RIP1  - [RIP1]·rc6 + [TNFR]·rc23
17  d  d t   P53  - [p53]·rc10 + [AKT]·rc18 - [p53]·rc27
18  d  d t   AKT  - [AKT]·rc18 + [HER2]·rc20
19  d  d t   STAT3  - [STAT3]·rc22
20  d  d t   TNFR  - [TNFR]·rc23
21  d  d t   DAPK  1  - [DAPK1]·rc28
22  d  d t   MOMP  +[DAPK1]·rc28 - [MOMP]·rc29 - [MOMP]·rc30 d  d t   A K T = - A K T   ⋅ r c  18   + H E R 2   ⋅ r c  20   ,   d  d t   p 53 = - p 53   ⋅ r c  10   + A K T   ⋅ | r c  18   - p 53   ⋅ r c  27   ,
where rc18 is the binding rate constant from AKT to p53 and [AKT] is the concentration of the AKT molecule. Also, rc20 represents the binding rate constant for HER2 to AKT. Similarly, for p53, variables in the corresponding equation incorporate information for the molecular concentrations and the binding rate constants for AKT and p53.
In our study, the ODEs were solved using the standard Matlab library function ode45. This function can interpret the differential equations in the form of rate reactions. The ode45 accepts the initial values of the molecules and the time duration in the form of function parameters [20,21]. Using this procedure, we derived the rate equations and time course data for all twenty-two signalling proteins. The initial concentrations of the signalling molecules were set to dimensionless numeric value of one [22]. Some of the rate constants were taken directly from the literature. These include Casp9 & Casp3, Casp9 & XIAP, Casp8 & BID, Casp3 & XIAP, Casp3 &Casp8, Casp3 & Apop, p53 & PUMA, Casp3 & Casp6 and XIAP & SMAC [4,23,24]. The remaining parameters were generated randomly but we verified them by Partial Least Square Regression (PLSR) test. We also performed sensitivity analysis of the model parameters by varying the model parameters to ensure overall robustness of the model to the small variations of parameters (see result section for PLSR test and sensitivity analysis).