<jupyter_start><jupyter_text>Cytoscape v3.5 CyREST Challenge Dynamic visualization of kinetic models in CytoscapeThis notebook contains a python workflow for the visualization kinetic model simulations within Cytoscape networks.In this workflow * kinetic models encoded in the Systems Biology Markup Language [SBML](http://www.sbml.org) are loaded* kinetic model simulations are performed * timecourse results are mapped on the corresponding Cytoscape networks which are read using the Cytoscaoe app `cy3sbml`. This workflow demonstrates how simulation results from computational models can be easily mapped on their model graph. By using `cy3sbml` simulation results in the graph are annotated with the information available from the SBML file, e.g. ontology information or metadata on parameters, species and reactions.The workflow uses a classical example model, the Elowitz Repressilator, with simulations resulting in an oscillatory change of the main players, three mRNA's and their respective protein levels.This workflow is valuable because it allows the automatic visualization of simulation results within the underlying biological network of the computational model. A wide range of models is available in SBML and the workflow can consequently be adapted easily to other use cases.* **solver**: <NAME> * **repository**: https://github.com/cytoscape/cytoscape-automation/tree/master/for-scripters/Python/Konig_SBML_Time_Course_-_2017_CyREST_ChallengeInformation on how to setup the python environment for executing this notebook is provided in:  https://github.com/cytoscape/cytoscape-automation/blob/master/for-scripters/Python/Konig_SBML_Time_Course_-_2017_CyREST_Challenge/README.md  Start Cytoscape 3 with cyrest & cy3sbml appsBefore executing this notebook start `Cytoscape>=3.5.1` with the following apps installed* `cyrest-app>=3.5.0` (http://apps.cytoscape.org/apps/cyrest)* `cy3sbml>=0.2.5` (http://apps.cytoscape.org/apps/cy3sbml)<jupyter_code>from __future__ import print_function, division
import pandas as pd
from matplotlib import pylab as plt
import roadrunner  
%matplotlib inline

# some settings for the plots
plt.rcParams.update({
        'axes.labelsize': 'large',
        'axes.labelweight': 'bold',
        'axes.titlesize': 'large',
        'axes.titleweight': 'bold',
        'legend.fontsize': 'small',
        'xtick.labelsize': 'large',
        'ytick.labelsize': 'large',
    })<jupyter_output><empty_output><jupyter_text>Load SBML modelIn the following we load the **Elowitz2000 model, a synthetic oscillatory network of transcriptional regulators** and user roadrunner, a high performance SBML simulator, to simulate the model. The model was downloaded from Biomodels as [BIOMD0000000012](https://www.ebi.ac.uk/biomodels-main/BIOMD0000000012).This model describes the deterministic version of the repressilator system. The authors of this model (see reference) use three transcriptional repressor systems that are not part of any natural biological clock to build an oscillating network that they called the repressilator. The model system was induced in Escherichia coli.In this system, **LacI (variable X is the mRNA, variable PX is the protein)** inhibits the tetracycline-resistance transposon **tetR (Y, PY describe mRNA and protein)**. Protein tetR inhibits the **gene Cl from phage Lambda (Z, PZ: mRNA, protein)**,and protein Cl inhibits lacI expression. **With the appropriate parameter values this system oscillates.**This model is described in the article:  A synthetic oscillatory network of transcriptional regulators.  <NAME>, <NAME>.  Nature. 2000 Jan; 403(6767):335-338  Abstract:Networks of interacting biomolecules carry out many essential functions in living cells, but the 'design principles' underlying the functioning of such intracellular networks remain poorly understood, despite intensive efforts including quantitative analysis of relatively simple systems. Here we present a complementary approach to this problem: the design and construction of a synthetic network to implement a particular function. We used three transcriptional repressor systems that are not part of any natural biological clock to build an oscillating network, termed the repressilator, in Escherichia coli. The network periodically induces the synthesis of green fluorescent protein as a readout of its state in individual cells. The resulting oscillations, with typical periods of hours, are slower than the cell-division cycle, so the state of the oscillator has to be transmitted from generation to generation. This artificial clock displays noisy behaviour, possibly because of stochastic fluctuations of its components. Such 'rational network design may lead both to the engineering of new cellular behaviours and to an improved understanding of naturally occurring networks.<jupyter_code># high performance SBML simulator
import roadrunner
sbml_model = "./models/BIOMD0000000012.xml"
r = roadrunner.RoadRunner(sbml_model)
model = r.getModel()<jupyter_output><empty_output><jupyter_text>Simulate SBML modelIn the following we perform a simple timecourse simulation which we than visualize in the `cy3sbml` Cytoscape network. The numerical results are stored in a `pandas DataFrame`. In the following the `sbml ids`, which we use for mapping on the Cytoscape network, and the first rows of the numerical solution are depicted.After this simulation step we have numerical time course data, which is used in the following to map on the network.<jupyter_code># reset the model
r.reset()

# select ids in numerical simulation
r.timeCourseSelections = ['time'] + model.getFloatingSpeciesIds() + model.getReactionIds()
# integration
data = r.simulate(start=0, end=200, steps=100)

# sbml ids from timeCourseSelections
sbml_ids = [name.replace("[", "").replace("]", "") 
            for name in r.timeCourseSelections]
print(sbml_ids)

# create pandas dataframe
df = pd.DataFrame(data[:, 1:], index=data[:,0], columns=sbml_ids[1:])
df.to_csv("./results/simulation.tsv", sep='\t')
df.head(10)<jupyter_output>['time', 'PX', 'PY', 'PZ', 'X', 'Y', 'Z', 'Reaction1', 'Reaction2', 'Reaction3', 'Reaction4', 'Reaction5', 'Reaction6', 'Reaction7', 'Reaction8', 'Reaction9', 'Reaction10', 'Reaction11', 'Reaction12']<jupyter_text>Plot numerical solutionNow we plot the solution to show the oscillating behavior of the system. mRNA and Protein levels of the system oscillate over time, with mRNA concentrations much lower than the protein amount (which can be seen on the log axis).<jupyter_code># The system oscillates
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(14, 7))
axes = (ax1, ax2)

for ax in axes:
    for sid in ['PX', 'PY', 'PZ', 'X', 'Y', 'Z']:
        ax.plot(df.index, df[sid], label=sid)
ax2.set_yscale('log')
for ax in axes:
    ax.set_ylabel('Value')
    ax.set_xlabel('Time')
    ax.legend()

plt.show()
fig.savefig("./results/simulation.png", bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Load SBML model in CytoscapeOn basis of the numerical results we do not understand the connection between the different mRNA and Proteins in this negative feedback model. For a better understanding what is going on we visualize the graph structure of the model.<jupyter_code>import json
import py2cytoscape
from py2cytoscape.data.cyrest_client import CyRestClient
from py2cytoscape.data.cynetwork import CyNetwork
from py2cytoscape.data.style import StyleUtil

# new cyrest client
cy = CyRestClient()

# reset Cytoscape session
cy.session.delete()
# load SBML network (this uses cy3sbml)
networks = cy.network.create_from('./models/BIOMD0000000012.xml')
print(networks)<jupyter_output>[<py2cytoscape.data.cynetwork.CyNetwork object at 0x7f6520059748>, <py2cytoscape.data.cynetwork.CyNetwork object at 0x7f651b32e6a0>, <py2cytoscape.data.cynetwork.CyNetwork object at 0x7f651b32e390>]<jupyter_text>After loading SBML models with cy3sbml three different networks are generated* All: contains all SBML information of the model (i.e. kinetics graph with units)* Kinetics: contains SBML information about the kinetics and relations to model components* Base: basic species-reaction graph of SBML (with additional GeneProducts if in the model)In the following we are mainly interested in the `species-reaction` graph of SBML which corresponds to * species: mRNA and protein* reactions: translation, transcription and mRNA/Protein degradationIn the following we select the base network<jupyter_code># Get information about SBML networks
base_net = None
for net in networks:
    # unique id of cytoscape objects
    net_suid = net.get_id()
    net_name = net.get_network_value(column='name')
    print(net.get_network_value(column='name'))
    if (net_name.startswith(b'Base')):
        base_net = net

# Selected the Base network (reaction - species graph of sbml file)
print('-'*60)
print (base_net)
print('-'*60)<jupyter_output>b'All: Elowitz2000 - Repressilator'
b'Base: Elowitz2000 - Repressilator'
b'Kinetic: Elowitz2000 - Repressilator'
------------------------------------------------------------
<py2cytoscape.data.cynetwork.CyNetwork object at 0x7f651b32e6a0>
------------------------------------------------------------<jupyter_text>Display the reaction-species graphHere we display the core model graph. With this figure it becomes more clear what is happening in the model and we clearly see the negative feedback cycle of the proteins on the transcription of the following proteins in the cycle.<jupyter_code># Adjust the node positions in Cytoscape and create the figure
from IPython.display import Image
Image(base_net.get_png())<jupyter_output><empty_output><jupyter_text>Mapping of simulation results to Cytoscape VisualStyle cy3sbml node tableAll cy3sbml table attributes are easily accessible from python. In the following we load the node table to find the nodes which correspond the the `sbml ids` in the simulation.<jupyter_code># cell = base_net.get_node_value(id=node0, column='name')
node_table = base_net.get_node_table()
node_table<jupyter_output><empty_output><jupyter_text>Find nodes for sbml idsWe now find the nodes which correspond to the species we simulated in the model.The suid (unique node identifiers) are retrieved and stored for lookup.<jupyter_code>node_suids = []

# logical indexing to find the nodes
for sid in sbml_ids[1:]:
    suid_index = node_table[node_table["sbml id"]==sid].index
    node_suids.append(suid_index.get_values()[0])

suid2sid = dict(zip(node_suids, sbml_ids[1:]))
suid2sid<jupyter_output><empty_output><jupyter_text>Map kinetic data to Visual StyleNow we map the kinetic simulation results on the VisualStyle of the Cytoscape network.We use the `NODE_HEIGHT` and `NODE_WIDTH` for mapping of the protein and mRNA concentrations.<jupyter_code># get network view
view_id_list = base_net.get_views() 
base_view = base_net.get_view(view_id_list[0], format='view')
print(base_view)
# create cytoscape images of timepoints
# the species amounts and reaction rates are mapped to the node sizes
from py2cytoscape.data.util_network import NetworkUtil as util
import time
import os
import numpy as np

# DataFrame for node views
df_vs_node = pd.DataFrame(index=node_suids, 
                  columns=['NODE_WIDTH', 'NODE_HEIGHT'])
# create cytoscape figure for timepoint
cypng_files = []
for timepoint in df.index:
    
    # set all node heights & widths
    for index, row in df_vs_node.iterrows():
        sbml_id = suid2sid.get(index)
        
        # Set node size from simulation results
        # area of nodes should be proportional to concentration, i.e. A^2 ~ H*W = H^2 = c
        row['NODE_WIDTH'] = np.sqrt(df.loc[timepoint, sbml_id])*3
        row['NODE_HEIGHT'] = np.sqrt(df.loc[timepoint, sbml_id])*3
        
        # normalisation with max value
        # row['NODE_WIDTH'] = 80 * df.loc[timepoint, sbml_id]/max(df[sbml_id])
        # row['NODE_HEIGHT'] = 80 * df.loc[timepoint, sbml_id]/max(df[sbml_id])

    # Apply style for timepoint
    base_view.batch_update_node_views(df_vs_node)
    
    cypng_data = base_net.get_png()
    cypng_file = "./results/images/cytime_{}.png".format(int(timepoint))
    cypng_files.append(cypng_file)
    
    if True:
        # write images
        with open(cypng_file, "wb") as f:
            f.write(cypng_data)
    else:
        # wait
        time.sleep(0.01)
# create kinetic plots for the respective timepoints
mplpng_files = []
for timepoint in df.index:
         
    # create matplotlib figure for timepoint
    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))

    for sid in ['PX', 'PY', 'PZ', 'X', 'Y', 'Z']:
        ax.plot(df.index, df[sid], label=sid, linewidth=2)
    ax.plot([timepoint, timepoint], [0,2500], label="__nolabel__", linewidth=2, color="black")
    ax.set_yscale('log')
    ax.set_ylabel('Value')
    ax.set_xlabel('Time')
    ax.legend()
    mplpng_file = "./results/images/mpltime_{}.png".format(int(timepoint)) 
    
    fig.savefig(mplpng_file, bbox_inches='tight')
    plt.close(fig)
    mplpng_files.append(mplpng_file)
# Combine kinetic plots and cytoscape images
from PIL import Image as PImage
compng_files =[]
for timepoint in df.index:
    # put smaller image in the center of cytoscape circle
    mplpng_file = "./results/images/mpltime_{}.png".format(int(timepoint)) 
    cypng_file = "./results/images/cytime_{}.png".format(int(timepoint))
    compng_file = "./results/images/comtime_{}.png".format(int(timepoint))

    background = PImage.open(cypng_file)
    foreground = PImage.open(mplpng_file)
    
    # calculate the inline position
    pos_x = int(1.0*background.width/2 - 1.0*foreground.width/2)
    pos_y = int(1.0*background.height/2 - 1.0*foreground.height/2)

    background.paste(foreground, (pos_x, pos_y), foreground)
    # background.show()
    background.save(compng_file)
    compng_files.append(compng_file)
# Create gif
import imageio
images = []
for filename in compng_files:
    images.append(imageio.imread(filename))
imageio.mimsave('./results/cytoscape-repressilator.gif', images)<jupyter_output><empty_output><jupyter_text>Display timecourseNow we have finished the dynamic visualization of the time course in Cytoscape. The following gif shows the resulting model behavior.   Display example imageIn the following is one of the example timepoints shown. The three proteins oscillate and peak in turns. As long as translation>degration for proteins the amount increases, when translation<degradation the amount of protein decreases.The mRNA oscillates also but with much lower amounts.<jupyter_code>from IPython.display import Image
Image('./results/images/comtime_100.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09449409838302518, 'fraction_numerical': 0.016237975028996383, 'mean_word_length': 4.798259493670886, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 53, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '54312', 'license': 'CC0-1.0', 'n_tokens_mistral': 4576, 'n_tokens_neox': 4178, 'n_words': 1513, 'path': 'for-scripters/challenge-2017/Konig_SBML_Time_Course_Python/cyrest-challenge.ipynb', 'repo_name': 'kozo2/cytoscape-automation', 'size': '14672'}