{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simulation visualization\n",
"\n",
"<div class=\"admonition note\">\n",
"**Topics**: Plotting, Mesh visualization\n",
"</div>\n",
"\n",
"The corresponding python script: [STEPS_Tutorial_Visual.py](https://github.com/CNS-OIST/STEPS_Example/tree/master/user_manual/source/API_2/scripts/STEPS_Tutorial_Visual.py)\n",
"\n",
"In a [previous chapter](STEPS_Tutorial_DataSaving.ipynb), we introduced how to save data to HDF5 files and generate XDMF files in order to visualize the data in the [Paraview](https://www.paraview.org/) software.\n",
"\n",
"In this chapter, we will use the STEPS visualization toolkit to interactively visualize STEPS simulations. To demonstrate this, we will use the simple IP3 receptor model introduced in the [surface-volume reactions chapter](STEPS_Tutorial_IP3.ipynb). The STEPS visualization toolkit was originally described in the following article: [Python-based geometry preparation and simulation visualization toolkits for STEPS](http://journal.frontiersin.org/Journal/10.3389/fninf.2014.00037/abstract).\n",
"\n",
"A full decription of all available classes and functions is available in the [API reference page](API_visual.rst).\n",
"\n",
"## Prerequisites\n",
"\n",
"The following third-party python packages are required to use the visualization toolkit.\n",
"\n",
"1. [NumPy](http://www.numpy.org/)\n",
"2. [PyQt](http://www.riverbankcomputing.co.uk/software/pyqt/download)\n",
"3. [PyQtGraph](http://www.pyqtgraph.org/)\n",
"\n",
"## IP3 receptor model visualization\n",
"\n",
"### Biochemical model\n",
"\n",
"We will use a spatial extension of the IP3 receptor model that we used in a [previous chapter](STEPS_Tutorial_IP3.ipynb). The main difference consists in adding diffusion rules for `Ca` and `IP3`:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import steps.interface\n",
"\n",
"from steps.model import *\n",
"from steps.geom import *\n",
"from steps.rng import *\n",
"from steps.sim import *\n",
"from steps.saving import *\n",
"from steps.visual import *\n",
"\n",
"# Source:\n",
"# Allbritton, N.L., Meyer, T., and Stryer, L. (1992). \n",
"# Range of messenger action of calcium ion and inositol \n",
"# 1,4,5-triphosphate. Science 258, 1812-1815.\n",
"DCST_Ca = 0.065e-9\n",
"DCST_IP3 = 0.283e-9\n",
"\n",
"mdl = Model()\n",
"r = ReactionManager()\n",
"\n",
"with mdl:\n",
" # chemical species objects\n",
" Ca, IP3, R, RIP3, Ropen, RCa, R2Ca, R3Ca, R4Ca = Species.Create()\n",
"\n",
" ssys = SurfaceSystem.Create()\n",
" with ssys:\n",
" # IP3 and activating Ca binding\n",
" R.s + IP3.o <r[1]> RIP3.s\n",
" RIP3.s + Ca.o <r[2]> Ropen.s\n",
" r[1].K = 1000000000.0, 25800.0\n",
" r[2].K = 8000000000.0, 2000.0\n",
"\n",
" # Inactivating Ca binding\n",
" R.s + Ca.o <r[3]> RCa.s\n",
" RCa.s + Ca.o <r[4]> R2Ca.s\n",
" R2Ca.s + Ca.o <r[5]> R3Ca.s\n",
" R3Ca.s + Ca.o <r[6]> R4Ca.s\n",
" r[3].K = 8889000.0, 5.0\n",
" r[4].K = 20000000.0, 10.0\n",
" r[5].K = 40000000.0, 15.0\n",
" r[6].K = 60000000.0, 20.0\n",
"\n",
" # Ca ions passing through open IP3R channel\n",
" Ca.i + Ropen.s <r[1]> Ropen.s + Ca.o\n",
" # Corresponds to Ca input ~ 20000/ms for open receptor\n",
" r[1].K = 8000000.0, 8000000.0\n",
" \n",
" vsys = VolumeSystem.Create()\n",
" with vsys:\n",
" # Create diffusion rules\n",
" Ca_diff = Diffusion.Create(Ca, DCST_Ca)\n",
" IP3_diff = Diffusion.Create(IP3, DCST_IP3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then load the tetrahedral mesh from file. This mesh already contains the patch and the compartments that were declared in the original IP3R model. We will then create the corresponding simulation object, and initialize the simulation:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Model checking:\n",
"No errors were found\n"
]
}
],
"source": [
"mesh = TetMesh.Load('meshes/ip3r_mesh')\n",
"\n",
"# Create random number generator\n",
"rng = RNG('mt19937', 512, 456)\n",
"\n",
"# Create reaction-diffusion solver object\n",
"sim = Simulation('Tetexact', mdl, mesh, rng)\n",
"\n",
"sim.newRun()\n",
"\n",
"# Setup initial condition\n",
"sim.cyt.Ca.Count = 1\n",
"sim.cyt.IP3.Conc = 2.5e-06\n",
"sim.ER.Ca.Conc = 0.00015\n",
"sim.memb.R.Count = 16\n",
"\n",
"rs = ResultSelector(sim)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In contrast with other chapters, we declare a `ResultSelector` object but we will not save any data with it. It will instead be used to describe which data we want to plot.\n",
"\n",
"### Declaring visualization elements\n",
"\n",
"The next step is to describe what we want to display:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Create control\n",
"sc = SimControl(end_time = 1.0, upd_interval = 0.0001)\n",
"\n",
"with sc:\n",
" # Plots\n",
" plots_d = PlotDisplay('IP3 Receptor Model')\n",
" with plots_d:\n",
" \n",
" TimePlot(rs.cyt.Ca.Conc, \n",
" title='Ca_cyt', \n",
" pen=(255, 0.647 * 255, 0), \n",
" data_size=1000, \n",
" y_range=[0, 15e-6], \n",
" y_label=('Concentration', 'M')\n",
" )\n",
" \n",
" TimePlot(rs.memb.Ropen.Count,\n",
" title='Ropen_memb', \n",
" pen=(255, 0, 255), \n",
" data_size=1000, \n",
" y_range=[0, 10], \n",
" y_label=('Count', '#')\n",
" )\n",
" \n",
" NewRow()\n",
" \n",
" SpatialPlot(rs.TETS(mesh.cyt.tets).Ca.Count,\n",
" title='Ca_cyt distribution anlong y-axis', \n",
" axis=[0, 1, 0], \n",
" y_range=[0, 30]\n",
" )\n",
" \n",
" TimePlot(rs.memb.MATCH('R.*').Count,\n",
" title='Receptor states', \n",
" data_size=1000, \n",
" y_range=[0, 17], \n",
" y_label=('Count', '#')\n",
" )\n",
" \n",
" # 3D displays\n",
" ER_d, CytIP3_d, CytCa_d, memb_d, full_d = SimDisplay.Create(\n",
" 'ER', 'Cyt IP3', 'Cyt Calcium', 'memb', 'Full view'\n",
" )\n",
"\n",
" with ER_d:\n",
" # Static mesh element\n",
" ElementDisplay(rs.ER, color=[0.678, 1, 0.184, 0.05])\n",
" # Dynamic element\n",
" ElementDisplay(rs.ER.Ca, color=[1, 0.647, 0, 1], spec_size = 0.005)\n",
"\n",
" with CytIP3_d:\n",
" ElementDisplay(rs.cyt, color=[0.941, 1, 0.941, 0.05])\n",
" ElementDisplay(rs.cyt.IP3, color=[1, 0, 0, 1], spec_size = 0.005)\n",
"\n",
" with CytCa_d:\n",
" ElementDisplay(rs.cyt, color=[0.941, 1, 0.941, 0.05])\n",
" ElementDisplay(rs.cyt.Ca, color=[1, 0.647, 0, 1], spec_size = 0.005)\n",
"\n",
" with memb_d:\n",
" ElementDisplay(rs.memb, color=[1, 0.973, 0.863, 0.05])\n",
"\n",
" # Different colors depending on receptor state\n",
" spec2Col = {\"R\" : [0.0, 0.0, 1.0, 1.0], \"RIP3\" : [1.0, 0.0, 1.0, 0.2], \n",
" \"Ropen\" : [1.0, 0.0, 1.0, 1.0], \"RCa\" : [0.0, 0.0, 1.0, 0.8], \n",
" \"R2Ca\" : [0.0, 0.0, 1.0, 0.6], \"R3Ca\" : [0.0, 0.0, 1.0, 0.4], \n",
" \"R4Ca\" : [0.0, 0.0, 1.0, 0.2]\n",
" }\n",
" ElementDisplay(\n",
" rs.memb.MATCH('R.*'), spec_size = 0.01,\n",
" color=lambda spec: spec2Col[spec.name]\n",
" )\n",
"\n",
" # Merge all sub-displays to the full one\n",
" full_d.merge(ER_d, CytIP3_d, CytCa_d, memb_d)\n",
" \n",
"# Enter visualization loop\n",
"sc.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to visualize our simulation, we create a `SimControl` object which represents the simulation control window and will be used to group all other windows. These latter need to be declared inside the `with sc:` block, in the same way as compartments are declared in a `with mesh:` block.\n",
"\n",
"After having declared all display windows, we need to launch the simulation with `sc.run()`. It is worth noting that the simulation object itself should not be controlled explicitely, it will be controlled by the simulation control window.\n",
"\n",
"There are two types of display windows that can be used:\n",
"\n",
"- **Plot displays** that contain plots showing the time evolution or spatial distribution of chosen quantities ; they are created with the `PlotDisplay` class\n",
"- **Simulation displays** that show 3D representations of the tetrahedral mesh with point representation of species ; they are created with the `SimDisplay` class\n",
"\n",
"In this guide, we will introduce how to create and populate these windows but we will not describe in details all their parameters, more information can be found in the corresponding [API reference page](API_visual.rst).\n",
"\n",
"### Plotting data\n",
"\n",
"A plot display window is created by constructing a `PlotDisplay` object, it takes a string corresponding to the window title as first parameter. Sub-elements in the window are then declared inside a `with` block that uses the display window:\n",
"\n",
"```python\n",
"with sc:\n",
" # PlotDisplays need to be declared in a SimControl\n",
" plotWin = PlotDisplay('My title')\n",
" with plotWin:\n",
" # Declare plots\n",
" # ...\n",
"```\n",
"\n",
"There are then two types of plots that we can add to the plot display:\n",
"\n",
"- **Time plots** showing the time evolution of chosen quantities ; they are created with the `TimePlot` class\n",
"- **Spatial plots** that can be used to display concentration gradients across space ; they are created with the `SpatialPlot` class\n",
"\n",
"#### Time plots\n",
"\n",
"In order to declare time plots, we need to provide at least one parameter to the `TimePlot` constructor: a result selector path to the data that should be plotted. The following code plots the same data as our main example, but without any of the optional formatting parameters:\n",
"\n",
"```python\n",
"with sc:\n",
" # ...\n",
" with plots_d:\n",
" TimePlot(rs.cyt.Ca.Conc)\n",
" TimePlot(rs.memb.Ropen.Count)\n",
" # ...\n",
" TimePlot(rs.memb.MATCH('R.*').Count)\n",
" # ...\n",
"```\n",
"\n",
"Note that, in contrast with conventional data saving, we do not need to add these result selectors to the simulation with the `toSave` method. It is possible to plot several lines in the sane `TimePlot` by passing a result selector path that encompasses several values. This is the case for the 3rd `TimePlot` in which we want to plot the count of IP3 receptors for each of its possible states.\n",
"\n",
"In the main example, several formatting keyword arguments are used to set the plot title (`title=...`), set the color of the lines (`pen=...`), limit the amount of data that can be displayed at a time (`data_size=...`), set the range of the y-axis (`y_range=...`), and set the label of the y-axis (`y_label=...`). Additional formatting parameters are described in the [API reference page](API_visual.rst#steps.API_2.visual.TimePlot).\n",
"\n",
"#### Spatial plots\n",
"\n",
"Spatial plots can be used to display the spatial distribution of some species along an axis. The following code corresponds to the same data as the main example, but without any of the optional parameters:\n",
"\n",
"```python\n",
"with sc:\n",
" # ...\n",
" with plots_d:\n",
" SpatialPlot(rs.TETS(mesh.cyt.tets).Ca.Count, axis=[0, 1, 0])\n",
"```\n",
"\n",
"The first parameter of the `SpatialPlot` constructor needs to be a result selector built using `TETS(...)` for tetrahedrons or `TRIS(...)` for triangles. For each of the tetrahedrons, it will query the corresponding value (number of 'Ca' here) and bin them according to the spatial location of the tetrahedron. The barycenter of the tetrahedrons is projected on the `axis` vector given as parameter. In our example, we project the locations on the `[0, 1, 0]` vector, i.e. the y axis. The locations are then binned into a number of bins that can be set with the `nbins` keyword parameter; its default value is 20.\n",
"\n",
"The spatial plot can then operate in two different modes, that can be selected with the `mode` keyword parameter:\n",
"\n",
"- **'distr'**, the default mode, for which the values associated with each tetrahedrons are being summed in each bin. When used with a `.Count` result selector, the resulting plot shows the distribution of species projected on the axis.\n",
"- **'mean'**, for which the values associated with each tetrahedrons are being averaged in each bin. When used with a `.V` result selector, the resulting plot shows the average electrical potential projected on the axis.\n",
"\n",
"Additional formatting parameters are described in the [API reference page](API_visual.rst#steps.API_2.visual.SpatialPlot).\n",
"\n",
"#### Layout\n",
"\n",
"Plots are added on a row from left to right, when `NewRow()` is called, a new row is created and subsequent plots will be added to this new row.\n",
"\n",
"### Visualizing mesh elements\n",
"\n",
"Mesh elements are displayed in simulation display windows which are created by constructing a `SimDisplay` object. Mesh elements are then attached to the window by declaring them inside a `with` block that uses the simulation window:\n",
"\n",
"```python\n",
"with sc:\n",
" # ...\n",
" simWin = SimDisplay.Create('My title')\n",
" with simWin:\n",
" ElementDisplay(rs.compA)\n",
" ElementDisplay(rs.compA.specB)\n",
" # ...\n",
"```\n",
"\n",
"In this example, we create a simulation window that will display the `compA` compartment as well as all `specB` species in this compartment. The surface of the compartment will be displayed in the window and the species will be represented as points. Note that compartments, patches, or simulation objects (species, channels, etc.) are all added to a simulation window with the same `ElementDisplay` class.\n",
"\n",
"#### Setting the color of elements\n",
"\n",
"In our main example, we use the `color=...` keyword parameter to set the colors of both compartments and species. Note that it takes a 4 elements list `[r, g, b, a]` with `r`, `g`, and `b` the colors itself (with each component being between 0 and 1) and `a` the alpha / transparency level, 0 meaning fully transparent and 1 fully opaque.\n",
"\n",
"The `color=...` keyword parameter can also take a function, as is the case in our main example for displaying the different receptor states. The lambda function we use there takes a specie as parameter and returns a color as a 4 elements list.\n",
"\n",
"The same thing can be done for setting the color of complexes according to their state. By default, if an element display is created with a path pointing to a complex, all states of the complex will be displayed in the same color. If the `color=...` keyword parameter is a function, it should take a complex state as parameter and return a color.\n",
"\n",
"#### Addtional formatting parameters\n",
"\n",
"The size of dots used to represent species can be adjusted with the `spec_size=...` keyword parameter. Additional keyword parameters are described in the [API reference](API_visual.rst#steps.API_2.visual.ElementDisplay).\n",
"\n",
"## Screenshots\n",
"\n",
"The following images correspond to the IP3 receptor example we presented in this chapter. They illustrate the different type of visualizations that are available in STEPS:\n",
"\n",
"### `PlotDisplay` - `TimePlot` and `SpatialPlot`\n",
"\n",
"<img src=\"images/visual_plots.png\"/>\n",
"\n",
"### `SimDisplay` - `ElementDisplay`\n",
"\n",
"#### Calcium in cytosol\n",
"\n",
"<img src=\"images/visual_Ca_Cyt.png\"/>\n",
"\n",
"#### Calcium in ER\n",
"\n",
"<img src=\"images/visual_Ca_ER.png\"/>\n",
"\n",
"#### IP3 in cytosol\n",
"\n",
"<img src=\"images/visual_IP3.png\"/>\n",
"\n",
"#### IP3 receptors on ER surface\n",
"\n",
"<img src=\"images/visual_IP3R.png\"/>\n",
"\n",
"#### Combined view\n",
"\n",
"<img src=\"images/visual_combined.png\"/>\n",
"\n",
"### `SimControl`\n",
"\n",
"<img src=\"images/visual_control.png\"/>"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}