The toy model presented here is a simple example consisting of three nodes: A, B and C. A can be activated by C and inhibited by B, B is only activated by A and C is activated by both A and B. This example illustrates the case of a simple loop in the transition graph. We propose to study two cases with two different transition rates for one of the transitions: fast and slow.

Since the model is small enough, we first chose to run it with asynchronous updating strategy in GINsim to visualize the transition graph.

The transition graph shows the fixed point [ABC]=[000]. It also shows a possible cycle: [101]=>[111]=>[011]=>[001]. At [001], the system can either continue to cycle or escape the cycle and go to the steady state [000]. We will explore the dynamics in the case that this "escape" rate is slow or fast.

The model can be written as follows:

Node A {
rate_up=(C AND (NOT B))  ? $Au : 0.0;
rate_down= B ? $Ad : 0.0;
}

Node B {
rate_up= A ? $Au : 0.0;
rate_down = A ? 0.0 : $Ad ;
}

Node C{
rate_up=0.0;
rate_down=((NOT A) AND (NOT B)) ? $escape : 0.0 ;
}

In the configuration file, we set the values for the parameters: Au, Bu, Ad, Bd and escape. The parameter "escape" corresponds to the degradation of C. In the fast case, escape=10, and in the slow case, escape=0.00001. Initial values for these simulations are set to 1.

BND and configuration files can be found here:

• BND file
• Fast escape rate Configuration file
• Slow escape rate Configuration file
• Slow escape rate Configuration file, long trajectory

Command lines:

  MaBoSS -c Four_cycle_FEscape.cfg -o Four_cycle_FEscape Four_cycle.bnd
  MaBoSS -c Four_cycle_SEscape.cfg -o Four_cycle_SEscape Four_cycle.bnd
  MaBoSS -c Four_cycle_SEscape_long.cfg -o Four_cycle_SEscape_long Four_cycle.bnd

There are several outputs of the simulations. Two of them present a particular interest for the biological interpretation of these simulations. The first one allows to follow the probabilities of each state and the entropy and transition entropy as a function of time (see probtraj files). The second one describes the asymptotic behaviour of the system, i.e. its possible stationary distributions (statdist files).

For the case of the fast escape, the probability to reach the stable steady state tends to 1. Both entropy and transition entropy tend to 0

For the case of the slow escape, if we only look at early response, the probability simulations do not show the presence of the stable steady state but rather the system seems stuck in a limit cycle. It is only after a long time that the system reaches the stable steady state (entropy and transition entropy reach 0). The two simulations show both the transient and the asymptotic behaviours.

Slow escape, early response:

Slow escape, late response:

The output files can be found below:

• Fast escape rate, probtraj
• Fast escape rate, statdist
• Slow escape rate, probtraj
• Slow escape rate, statdist
• Slow escape rate, probtraj, long trajectory
• Fast escape rate, statdist, long trajectory
• Excel file

The model is based on a published Boolean model of mammalian cell cycle (PMID:16873462), available on GINsim repository. It describes the signaling network that controls the mammalian cell cycle.

The model is described in the bnd file below. In the description of the model, we have included the possibility to simulate mutants. In the rate_down and rate_up terms, the parameters related to each mutant are referred to as $gene_del for deletion. They are set to 0 in the configuration file for the wild type and set to 1 when a mutant is modeled.

The corresponding BND and configuration files can be found below:

• BND file
• Configuration file, to be tuned for mutant models, initial condition corresponding to cell cycle start
• Configuration file, to be tuned for mutant models, random initial condition for studying stationary distribution

Command lines:

  MaBoSS -c cellcycle_runcfg.cfg -o cellcycle_runcfg cellcycle.bnd
  MaBoSS -c cellcycle_runcfg_randinit.cfg -o cellcycle_runcfg_randinit cellcycle.bnd

Damped oscillation can be observed on probablities of active cyclin.

Modeling with a random inital condition produces two clusters, that represent two indecomposable stationary distributions

The first cluster can be interpreted in term of a disynchronized proliferating cell population

The second cluster can be interpreted as a fixed point, representing the G1 phase.

The corresponding output files can be found below:

• probtraj file, initial conditition corresponding to cell cycle start
• statdist file, random initial conditition
• excel session file for probability trajectories
• excel session file for stationary distribution

The cell fate decision example is based on a published Boolean model (PMID:20221256). It describes how the cell chooses between three different cell fates (apoptosis, NFkappaB activation and necroptosis) in response to cell death receptor engagement (through Fas and TNFR).

The model is described in the bnd file below. In the description of the model, we have included the possibility to simulate mutants. In the rate_down and rate_up terms, the parameters related to each mutant are referred to as $gene_del for deletion and $gene_oe for over-expression respectively. They multiply a constant (1000 here) and are set to 0 in the configuration file for the wild type and set to 1 when simulated.

The corresponding BND and configuration files can be found below:

• BND file, including mutants as external variables
• Configuration file, to be tuned for mutant models

Command line:

  MaBoSS -c CellFateModel.cfg -o CellFateOut_wt CellFateModel.bnd 

For the wild type case, there are three stable steady states distributed as follows:

All the trajectories form three clusters that can be compared to the results of the published model. Qualitatively, the reachable phenotypes are the same. These probabilities are not absolute values as they depend on the computation and the details of the model itself. In the published model, the trajectories were computed over the transition graph.

The published model showed the following phenotypes, where A stands for Apoptosis, N for Necrosis and S for survival:

The model simulated in MaBoSS shows the same phenotypes with slightly different distributions:

Where the clusters correspond to:

[Apoptosis]=[TNF--TNFR--FADD--DISC_TNF--CASP8--CASP3--BAX--ATP--MOMP--SMAC--Cyt_c--apoptosome--Apoptosis]
[Necrosis]=[TNF--TNFR--FADD--DISC_TNF--CASP8--BAX--MPT--ROS--MOMP--SMAC--Cyt_c--NonACD]
[Survival]=[TNF--TNFR--FADD--DISC_TNF--RIP1--cIAP--RIP1ub--RIP1K--IKK--NFkB--cFLIP--ATP--XIAP--Survival]
	  

The simulation of CASP8 deletion shows the disappearance of one phenotype: apoptosis (in accordance with published experiments):

The corresponding output files can be found below:

• probtraj of wild type
• statdist of wild type
• probtraj with CASP8 deletion
• statdist with CASP8 deletion
• probtraj with XIAP deletion
• statdist with XIAP deletion
• excel file with mutant outputs

We consider a published model of p53 response to DNA damage (PMID:19233211). p53 interacts with Mdm2, which appears in two forms, cytoplasmic and nuclear. On one hand, p53 upregulates the level of cytoplasmic Mdm2 which is then transported into the nucleus and inhibits the export of nuclear Mdm2. On the other hand, Mdm2 facilitates the degradation of p53 through ubiquitination. In the model, stress regulates the level of DNA damage, which in turn participates in the degradation process of Mdm2. p53 inhibits the DNA damage signal by promoting DNA repair.

The transition graph contains two cycles and a fixed point [p53 Mdm2C Mdm2N Dam] = [0010] where nuclear Mdm2 is on and the rest is off.

BND and configuration files:

• BND file
• Configuration file

Command line:

MaBoSS -c p53_Mdm2_runcfg.cfg -o p53_Mdm2 p53_Mdm2.bnd
	  

The trajectories of the probabilities of p53 being equal to 1 and p53 being equal to 2, with the initial condition: [p53 Mdm2C Mdm2N Dam] = [0*11] are plotted below:

And the entropies:

The corresponding output files can be found below:

• probtraj
• statdist
• excel file

The model of segment polarity genes patterning is based on a published model (PMID:12782112). The following analysis is based on another publication (PMID:20655356); in particular, this publication studied the effect of random noise applied to the signaling network of segment polarity genes. The model consists of four identical signaling pathways, that are connected to each other. It represents the last step of drosophila patterning: the expression of the segment polarity genes. More precisely, this last step consists of a finer spatial patterning, driven by a set of genes (the segment polarity genes), whose expression last for a long time (compare to other genes involved in Drosophila segmentation).

Like in the two published papers, the model consists of four signaling pathways, one for each adjacent cell. mRNA and proteins are considered separately. Cell-cell interactions are represented by connections of one cell to the adjacent ones. For example, gene EN is transcriptionally activated by WG through the membrane; this is translated by including the following connections: WG4->en1, WG2->en1, WG1->en2, WG3->en2, WG2->en3, WG4->en3, WG3->en4, WG1->en4. The node en is written in lower case, because it represents the mRNA associated to the protein EN. Note that the model is "periodic", because cell 1 is adjacent to cell 4. SMO and PH are replaced by direct connections.

The logical model is translated in MaboSS language. Two values of transition rates are considered: one for transcriptional influences ($rate_RNA=1), one for post-transcriptional influences ($rate_protein=10). The effect of random perturbation on every node has been included in the BND file; it can be tuned by the parameter $noise.

The expected biological patterning is represented by the reference state below; none of the four cells are identical in term of node states. As described in the publications above, the model should produce the expected patterning, starting from a given initial state, when the pair rule gene SLP has the following pattern: SLP1,2=0 and SLP3,4=1. Two cases are considered: no noise ($noise=0) and low noise ($noise=0.01).

The corresponding BND and configuration files can be found below:

• BND file
• Configuration file with no noise
• Configuration file with noise
• Configuration file with noise, tuned for computing stationary distribution

A specific configuration file for computing stationary distribution is provided for the following reason: it is easy to show that random noise on nodes expression ($noie>0) implies that there is only one indecompsable stationary distribution, with all network states with a non-zero probability (for the given intial condition). Therefore, the clustering algorithm for identifying indecomposable stationary distributions cannot give good results, seprating artificially stationary distribution estimates. In that case, it is more efficient to increase the maximum time, to lower the set of trajectories and to set the threshold for clustering to zero.

Command lines:

MaBoSS -c D_Pattern_ContT_noIntNode.cfg -o D_Pattern_ContT_noIntNode D_Pattern_ContT.bnd
MaBoSS -c D_Pattern_ContT_noIntNode_noise.cfg -o D_Pattern_ContT_noIntNode_noise D_Pattern_ContT.bnd
MaBoSS -c D_Pattern_ContT_noIntNode_noise_4statdist.cfg -o D_Pattern_ContT_noIntNode_noise_4statdist D_Pattern_ContT.bnd
	  

When there is no noise on nodes, the expected patterning is reached. This can be seen on the time trajectories of the Hamming distance distribution (with expected patterning as the reference state):

The asymptotic behavior consists of fixed points. Their respective probabilities are:

The blue fixed point corresponds to the correct patterning.

When there is noise, the expected patterning is no longer stable. This can be seen on the time trajectories of the Hamming distance distribution (with expected patterning as the reference state):

In that case, the asymptotic behavior consists of a single indecomposable stationary distribution:

The most probable state (in yellow) corresponds to a "broad type" patterning (this state is defined in the publications above): two adjacent cells have identical node states.

The corresponding output files can be found below:

probtraj file with no noise

statdist file with no noise

probtraj file with noise

statdist file with noise

This model was used as an example for the standard language SBML qual and published in 2013 by Chaouiya et al. The model can be found in BioModels database with the reference BIOMD0000000562.xls

$ MBSS_FormatTable.pl chaouiya_maboss.bnd chaouiya_maboss.cfg