# 3.1. Example combining SBMLqual, CNORdt, CNORode, CNORfeeder packages¶

This tutorial is for advanced users. It combines the **CellNOptR**, **CNORode**
and **CNORfeeder** packages. The goal of this tutorial is to provide the material to reproduce some results shown or discussed about in the SBMLqual paper, which is available in the arxiv.

The *first part* reproduces one of figure shown in the paper. This figure was reproduced by the different software used in the paper using the same
SBMLqual model as input. Here, we show how one can use CellNOptR oftware to simulate the same dynamical states as the one shown in the paper.

The *second part* is more about CellNOptR software itself: we deal with the creation of some data sets (with CNORode), optimisation (with CNORdt) and usage of CNORfeeder to find missing links.

Warning

In order to reproduce the examples and results provided here below, you must use some R packages available on cellnopt website rather than the current version available on BioConductor website, which are currently older than latest version required to reproduce the results. The examples were generated with these package/version:

- CellNOptR “1.8.2”
- CNORdt “1.5.1”
- CNORfeeder “1.2.1”
- CNORode “1.3.3”

## 3.1.1. Reading a SBMLqual and simulate using CNORdt¶

In CellNOptR, it is possible to load a SBMLqual file since version 1.7.8. Given the model, one can then decide to use different formalisms to simulate the dynamical states.

In this section, we provide a script called `2013/sbml/simulate.R`
that loads a SBMLqual and performs a simulation using another package called CNORdt. This package provide a discrete time formalism to simulate the dynamical states.

You can download the script from the link above and execute it as follows:

`R script --no-save --no-restore < simulate.R`

The SBMLqual file can be downloaded here:
`2013/sbml/data/ModelV5.xml`. The file is expected to be found in the
directory: **./data**. Otherwise you can change the script.

Each combination of stimuli lead to a different dynamical states. There are 4 combinations and therefore 4 pictures are generated. Here is one of them for which TNFa and EGF are stimulated:

## 3.1.2. Creating data from SBMLQual and optimising a PKN on the data set.¶

In this section, we will do the following:

- Using the same
`SBMLQual model`as above, which we will call**True model**, we first generate the dynamical states using the ODE formalism (CNORode). The states are saved in MIDAS format. - Perform an optimisation (using CNORdt formalism) of the True model against the data for sanity check (the optimised model should be identical to the True model !!) and they fit almost perfect.
- Then, we will create a Prior Knowledge Network (PKN) that is slightly different from the True model.
- Using the PKN model and the data, we optimise the PKN model using a CNORdt formalism. We will see that there are discrepancies due to the differences between the True model and the PKN, in particular that there is a missing link in the PKN.
- Finally, we will use CNORfeeder to infer the missing link based on the data. An optimisation will show that the CNORfeeder found the missing link.

This example is based on the PKN and True model available in http://iopscience.iop.org/1478-3975/9/4/045003

Here is the code to generate the image corresponding to the True model:

```
library(CellNOptR)
model = readSBMLQual("ModelV5.xml")
library(CNORdt)
data(CNOlistPB, package="CNORdt")
plotModel(model, CNOlistPB)
```

Note

the CNOlistPB variable is optional and used solely for setting colors on nodes

Note

In the network plot, you can see that (1) there are 2 AND gates; they do
not have any corresponding OR gates (2) there is a link from *traf2* to *mkk7*.

### 3.1.2.1. Creating the data¶

To create the data using an ODE formalism, you will need a set of parameters
available in `2013/sbml/data/paramsODE.RData`.

Warning

if the SBML is changed, the order of the species may be different an you would need to map the parameters stored in the RData file above with your new SBML file.

```
library(CNORode) # for the simulation
library(CNORdt) # for the data
data(CNOlistPB, package="CNORdt")
# Read the model with a dummy link on ph self loop
model = readSBMLQual("ModelV5.xml")
# load a variable called params that is required for the ODE run
load("paramsODE.RData")
timeSignals = seq(0,30,1)
CNOlistPB$timeSignals = seq(0,30,2) # a trick to keep the original dt
# somehow stimuli set to exactly zero leads to NA so, replace them with
# tiny value but must set back to zero later on
CNOlistPB$valueStimuli[1,1] = 0.01
CNOlistPB$valueStimuli[1,2] = 0.01
sim_data = plotLBodeModelSim(CNOlistPB, model, ode_parameters=params,
timeSignals=timeSignals, maxStepSize=0.035, show=F)
CNOlistPB$timeSignals = timeSignals
CNOlistPB$valueStimuli[1,1] = 0.
CNOlistPB$valueStimuli[1,2] = 0.
# finally converts the simdata to a cnolist to be saved in a file.
cnolist = simdata2cnolist(sim_data, CNOlistPB, model)
# keep only species of interest
species2select = c('raf1','erk', 'ap1', 'gsk3','p38', 'nfkb')
conditions2select = c(1,2,3,4,5,6,7,8,9,10)
for (i in seq_along(cnolist@signals)){
cnolist@signals[i][[1]] =
cnolist@signals[i][[1]][conditions2select,species2select]
}
plot(cnolist)
```

Finally, let us save the new data set into a file:

```
writeMIDAS(cnolist, "data.csv", overwrite=T)
```

### 3.1.2.2. Checking¶

Let us first check that the model optimised against the data gives back the SBMLqual model.

```
library(CNORdt)
model = readSBMLQual("ModelV5.xml")
cnolist = CNOlist("data.csv")
# test that we can retrieve the exact data/model with itself
upperB = 10
lowerB= 0.8
boolUpdates = 30
bString=rep(1,30)
# optimised with itself to check and get the best RMSE
# Here we can set the maxtime to only 20 seconds
opt1 = gaBinaryDT(CNOlist=cnolist, model=model,
initBstring=rep(1,length(model$reacID)),
boolUpdates=boolUpdates,maxTime=20, lowerB=lowerB, upperB=upperB,
stallGenMax=20)
```

We can check that the best model is the model itself (all bits are on) and that the RMSE is small:

```
> opt1$bString
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> opt1$bScore
0.00142731
```

Finally, the following plot shows that errors across conditions and readouts are all small (white) by using this code:

```
cutAndPlotResultsDT(CNOlist=cnolist, model=model, upperB=upperB,
lowerB=lowerB, bString=opt1$bString, boolUpdates=boolUpdates)
```

### 3.1.2.3. Defining a PKN model¶

We use the original True model stored in the SBMLqual (image on top) but alter it in the following ways:

- We remove the 2 AND gates since we don’t know a priori the logical gates
and add back the OR gates.

We add the link tnfr–>pi3k

We add the link pi3k–>rac–>map3k1

We remove the link traf2->ask1->mkk7

The PKN is in SIF format available here `2013/sbml/PKN_noask1.sif`

### 3.1.2.4. Optimising the PKN against data using CNORdt¶

First, we can check indeed that the True model can be used as a PKN to be optimised against the data generated. In theory, the final model should be identical to the True model.

```
pknmodel2 = readSIF("PKN_noask1.sif")
model = preprocessing(cnolist, pknmodel2, compression=F, expansion=T)
plotModel(model, cnolist)
# let us try this one where the additional links
# TNFR->PI3K and PI3K-RAC-Mapk31 are off as well as the unneeded AND and OR
gates.
# If you run the optimisation long enough, you'll get the optimised model to
# have the following bitstring:
optbs = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,0,0,0,1,0,0,0,1,0)
# It gives an RMSE = 0.031
opt1 = gaBinaryDT(CNOlist=cnolist, model=model,
initBstring=optbs,boolUpdates=boolUpdates,maxTime=1000, lowerB=lowerB,
upperB=upperB, stallGenMax=200, elitism=2, popSize=100, sizeFac=1e-4)
```

Finally, the following plot shows that errors across conditions and readouts are all small (white) by using this code:

```
cutAndPlotResultsDT(CNOlist=cnolist, model=model, upperB=upperB,
lowerB=lowerB, bString=optbs, boolUpdates=boolUpdates)
```

The fit is not good for the AP1 species. It means that

- The additional links that were added have been correctly removed (e.g., tnfr->pi3k)
- The OR gates added have been removed while the correct required AND gates (only) were kept
- The missing link on ASK1 prevents the AP1 to be fitted correctly when TNFa is on.

### 3.1.2.5. CNORfeeder to infer new links¶

We can fix the issue of the missing link by using the CNORfeeder package to infer missing links based on the data solely.:

```
#
library(CNORfeeder)
pknmodel = readSIF("PKN_noask1.sif")
model = preprocessing(cnolist, pknmodel, compression=F, expansion=T)
# find new links
BTable = makeBTables(CNOlist=cnolist, k=2, measErr=c(0.1,0))
modelIntegr = mapBTables2model(BTable=BTable, model=model, allInter=F, compressed=FALSE)
modelIntegr$reacID[modelIntegr$indexIntegr]
plotModel(model=modelIntegr, CNOlist=cnolist, indexIntegr=modelIntegr$indexIntegr)
```

Note

we set compressed=False when calling mapBTables2model function to be in agreement with the preprocessing step

You can see here below, in purple, the links that have been inferred.

There is a proposed link between TNFa and ap1. This is related to the missing link we are looking for.

Does it improve the optimisation ?

```
# We could run the optimisation, which takes a while. The following
# bitstring is the optimal one
optbs_feeder = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1)
# we can check this hypothesis by using it as the initial guess:
opt1 = gaBinaryDT(CNOlist=cnolist, model=modelIntegr, boolUpdates=boolUpdates,maxTime=1000, lowerB=lowerB, upperB=upperB,
stallGenMax=20, elitism=4, popSize=50, sizeFac=1e-4, initBstring=optbs_feeder)
# and check the final model
plotModel(modelIntegr, cnolist, bString=optbs_feeder)
```

One can check with the CNORdt optimisation that the new MSE is down to 0.00148 as compared to 0.031. The RMSE obtained is equivalent to the one obtained with hte True model (0.00142). The difference is due to the data being generated with a pure ODE formalism and the discrete time being used for the optimisation. We can check on the fitness plot that AP1 is now optimised properly:

```
cutAndPlotResultsDT(CNOlist=cnolist, model=modelIntegr, upperB=upperB,
lowerB=lowerB, bString=optbs_feeder, boolUpdates=boolUpdates)
```

So, the CNORfeeder analysis tells us that there is a missing link between TNFa stimulus and AP1 readouts. We know from the PKN that indeed there is a cross talk via ASK1 species. Without that knowledge, one can use different resources (e.g., intact) to figure out existing interactions between TNFa and ap1.