This tutorial requires BioNetGen 2.2.5, RuleBender, MATLAB2013a, CellOrganizer 2.1, Blender 2.70a, CellBlender 1.0 RC3, and MCell 3.2.1.
Note
The CellOrganizer section of this tutorial is not supported on the Windows operating system. If you wish to use the Import SBML functionality in MacOSX or Windows you need to use the Repository Cellblender version, Python 2.7 and LibSBML 5.10.
In part one, you will receive an introduction to working with spatial modeling. Using CellBlender 1.0 RC3 and MCell 3.2.1, you will create and run a simple biochemical system and geometries. Lastly you will learn to import the geometries from SBML-spatial and check that these meshes are manifold and watertight. You will be using these complex geometries later.
Using BioNetGen 2.2.5 (BNG), you will learn to model biochemical systems using rule based modeling. You will use powerful modeling paradigm to create and analyze ODE simulations of signal transduction. You will then learn how to export biochemistry from BNG to Systems Biology Markup Language (SBML) and import it into CellBlender to create a simulation using the complex geometries you worked on in the previous section.
Using CellOrganizer 2.1 you will learn to train generative models of cellular organization from fluorescence microscopy images. Next you will learn to sample the parameter space of your cellular models to create geometries specific to your interests. You will then synthesize in silico instances of cellular organization using these models. Lastly we will demonstrate how to import an SBML file into CellOrganizer to automatically determine what models are required for the biochemical system within the SBML file and synthesize instances from those models that are exported to SBML-spatial. Lastly we will return to your CellBlender simulations you started the in part two to analyze the results of these spatially realistic simulations and the impact of spatial organization on these simulations.
Blender 2.70a, CellBlender, MCell
Note
You may have already done this next step if you've followed along with earlier tutorials.
First download and install Blender 2.70a. Next download CellBlender 1.0 RC3. Open Blender. Click on File->User Preferences. In the "Add Ons" tab select "Install from file" and select the cellblender.zip you downloaded previously. If you then search the list of Add Ons for CellBlender you will now see a CellBlender add on that can be activated by clicking the box on the right hand side of the list. To save this as your default click "Save User Settings" before exiting this window. You should now see several "CellBlender-" sections in the options on the right panel particularly in the "Scene" tab.
Harris et al. present in [1] a model of receptor-mediated signaling coupled with nuclear transport and transcriptional gene regulation. The full model is shown in :num:`Fig. #cbngfull`
Throughtout this tutorial we will show the user how to implement this model using CellBlender and various external tools. During this first section, we will implement the reduced version shown in :num:`Fig. #cbngreduced` which includes ligand-receptor binding, ligand binding, ligand-receptor dimerization and TF-receptor binding.
Full Model Diagram
Reduced Model
Our basic system is composed of an extracellular matrix a cytoplasm and its membrane. For this example model we will use a basic geometry composed of an icosphere enveloped by a cube in CellBlender (:num:`Fig. #simple`) Name the sphere CP (for cytoplasm) and the bounding box EC (for extracellular matrix). Alternatively you can use the MDL import feature and use a premade geometry (simpleGeometry.mdl).
Simple Geometry
Add the following reaction network:
L + R <-> L_R kp_LR,km_LR L + L <-> L_L kp_LL,km_LL TF + TF <-> TF_TF kp_TF_TF, km_TF_TF R + TF <-> R_TF kp_R_TF,km_R_TF
And the following parameters:
NaV 6.022e8 kp_LR 0.1*NaV km_LR 1 kp_LL 0.1*NaV km_LL 1 kp_TF_TF 0.1*NaV km_TF_TF 1 kp_R_TF 0.1*NaV km_R_TF 0.1
With molecules:
L (3D) diffusion: 8.52e-7 R (2D) diffusion: 3.32e-7 TF (3D) diffusion:8.52e-7
And the following release sites:
Molecule: L Object: EC[ALL] - CP[ALL] Quantity type: Concentration/Density Quantity: 1000 / (NaV * 27) Molecule: R Object: CP[ALL] Quantity type: Concentration/Density Quantity: 200/ 24 Molecule: TF Object: CP[ALL] Quantity type: Concentration/Density Quantity: 200/ (NaV *8)
Under model initialization select 1000 iterations. Save your project and run the simulation!
The reaction network we have implemented is a highly minimalistic version when compared to the full reaction network. In this exercise we have attempted to show you the difficulty of simulating highly complex biological systems. We have developed a tool that we will present to you in the following sections detailing a procedure for defining complex reaction networks and using realistic geometries in a structured way.
BioNetGen 2.2.5, RuleBender 2.0, JavaSE 6.
RuleBender is a graphical user interface for BioNetGen. It includes several features like an integrated development, analysis and simulation environment that you will learn to use during this tutorial. You can get the latest version of RuleBender (which includes a copy of BioNetGen) from http://www.rulebender.org.
This section introduces the following concepts:
- Molecule: BioNetGen’s basic simulation unit. It is a structured object containing compartments that can be modified or bind to each other
- Species: A complex containing one or more molecules
- Rule: Rules are the generators of species in a BioNetGen model. Through the use of patterns and wildcards, a single rule can correspond to several reactions describing the same chemical process occurring under different biological contexts.
template.bngl, lr.bngl (within bngFiles.zip)
To teach the student the basics of how to set up a basic BioNetGen file with parameters, molecule and reaction definitions and observables.
Open the template file we provided in RuleBender. A BioNetGen file is comprised of a model definition followed by actions. To begin the model definition block we start with the line
begin model
The model definition is comprised of five required blocks and two optional
ones. The required blocks are parameters, molecule types, species,
observables, and reaction rules. The optional blocks, which will be
covered later, are compartments and functions.
The first step is to define the parameters we will use in the definition of our model. These parameters can take on any numerical value, although parameters used to define concentrations or rate constants should be non-negative. The syntax is the following:
begin parameters L0 1000 # Initial number of L molecules R0 20 # Initial number of R molecules kp1 1 # Bimolecular rate constant for L-R binding (1/# 1/s) km1 1 # Unimolecular rate constant for L-R unbinding (1/s) end parameters
Note that it is possible to define parameters using expressions involving previously defined parameters. This is convenient for documenting how unit conversions were performed.
The next step is to establish the set of molecules we will use in our system. In this case, it is sufficient to define a ligand molecule (L) and a receptor (R), each with a cognate binding site, r and l respectively:
begin molecule types L(r,d) R(l) end molecule types
The next step is to specify the species that are initially present in the
system along with their intial concentrations, which is done in the species
block:
begin species L(r,d) L0 R(l) R0 end species
Here, the two species initially present are free ligand, an L molecule with
its r and d sites unbound, and free receptor, an R molecules with its l site
unbound. Their initial concentrations are set to R0 and L0
respectively.
Next we define the outputs of the models, which are called observables.
Observables are defined as sums over the concentrations of species selected by
a particular specified pattern or set of patterns. Patterns can be thought of
as search terms that require the matched species to have a particular set of
properties. An example observables block is:
begin observables Molecules L_tot L() Molecules L_free L(r) Molecules L_bound L(r!1).R(l!1) end observables
The first observable, L_tot, computes the total number of L molecules in
the system. It does this by matching every occurence of an L molecule
regardless of the state of its r component. This illustrates the ‘don’t write
don’t care’ principle in BioNetGen, which means that what we don’t include in a
pattern doesn’t affect the match. Here, the pattern matches an L but doesn’t
add any requirements on the component states, so all L molecules are matched by
the pattern. The second observable, L_free, uses a pattern that lists a
component a single component, r, and by so doing requires that the matching L
molecules has an r component that is unbound. The third observable,
L_bound, counts the number of L molecules that are bound to an R molecule
using a more complicated pattern that illustrates the syntax used to specify
bonds. The ’!’ after a component is used to refer the binding state, and the
numerical index that follows is used to tag the endpoints of a particular bond.
Here, the bond with index 1 links the r component of L to the l component of R.
Each bond should have exactly two enpoints. Different indices are used to
refer to different bonds. The scope of the bond indices in BNG is always a
single pattern.
Rules are the central component of a BNG model, as they define what reactions can take place in the system. In this simple model, the reversible binding of L and R is specified as
begin reaction rules L(r) + R(l) <-> L(r!1).R(l!1) kp1, km1 end reaction rules
The left hand side of the rule (the part before the arrow) defines the
reactants and the right hand side defines the products. Reactants and products
are each selected by patterns that may match many different species in the
systems, and thus each rule may generate many reactions. The arrow for a rule
may be either unidirectional (->) or bidirectional (<->). The rate of
each generated reaction is determined by the rate constant expression
(unidirectional) or pair of rate constant expressions (bidirectional) that are
listed following the patterns. For the models we will use in this tutorials,
all reactions will follow elementary reaction kinetics meaning that the rate is
simply the product of a rate constant, given by the specified parameter, and
the reactant species concentrations. This is the only type of reaction rate law
that is currently allowed in a model that will be imported into MCell.
Specifying the reaction rules completes our description of the model, which we now terminate with the line
end model
The remaining lines of the BNGL (BioNetGen language) file specify actions
that are to performed on the model. If the model is to be simulated using
either ODEs or the SSA, the first action that needs to be performed is to
generate the reaction network starting from the defined species (see
species block above) and the reaction rules, by issuing the following
command:
generate_network({overwrite=>1})}
Here, we have used the overwrite option to ensure that any previously
generated network with the same base name will be overwritten (not necessary if
you are using RuleBender). More details about options that can be passed to
BioNetGen action commands can be found here.
The second action we want to perform is a simulation of the model using ODEs, which we specify using
simulate({method=>"ode",t_end=>5,n_steps=>120})
This generates a single trajectory by solving the ODEs defined by the rate
equations for the generated reaction network on the time interval [0,5],
sampled at 120 points (not including the starting point). :num:`Fig. #plot1`
shows this trajectory. We can see how Lig\_Bound initially increases until the
system reaches equilibrium between the forward and backward reactions.
Time Trajectory of a Ligand-Receptor Model
A second way to make use of components in BNG is through the use of internal states. A modeler can use these, for example, to represent post-translational modifications at a particular site or conformational states of a particular domain. For example, we can add a phosphorylation site to the receptor molecule in our previous example by modifying the molecule type declaration:
begin molecule types L(r) R(l,Y~0~P) end molecule types
where we have indicated on the second line that the molecule R contains an
additional component Y (representing a tyrosine residue) that can be in either
the state 0 (representing unphosphorylated) or the state P
(representing phosphorylated).
To model a phosphorylation reaction, we add the rule
begin parameters p1 1 d1 0.1 end parameters begin reaction rules ... R(l!+,Y~0) -> R(l!+,Y~P) p1 end reaction rules
In this rule component T transitions from state 0 to state P. We
are also introducing some new concepts with component R(l!+). Fist, the
!+ symbol is used to indicate that we are requiring component R(l) to
be bound to some molecule without specifying which. Separately, in this rule
the pattern R(l!+) is being used as context for the rule (a
precondition). In other words, for molecule R to become phosporylated it is
necessary for its component R(l) to be bound to something.
The phosporylation rule was unidirectional, so we also need to define a dephosphorylation reaction. We will assume that dephosphorylation doesn’t have a contextual requirement - here it doesn’t require binding of the ligand molecule, so we have
begin reaction rules ... R(Y~P) -> R(Y~0) d1 end reaction rules
In order to define an observable that tracks the number of phosporylated tyrosine residues we use the following pattern syntax.
begin observables Molecules R_Phospo_unbound R(Y~P) Molecules R_Phospo R(Y~P!?) end observables
The R(Y~P)!? syntax is used to indicate that we wish to count all
R molecules with a phoporylated Y component regardless of its bound state.
In contrast, the R(Y~P) is used to indicate that we strictly wish to
count the cases were Y is unbound. In this particular example it bears no
difference since we have not specified any rule that defines Y to also be a
binding site.
Compartments in BioNetGen (cBNGL) allow us to explicitly model the compartmental organization of a cell. Consider our example model we referred to in the previous section. :num:`Fig. #compartment` shows an abstract view of the hierchical structure of the model. cBNGL syntax allows us to directly use this hierarchy inside our model as such:
begin parameters vol_EC 20.0 vol_CP 4.0 vol_EN 0.5 sa_PM 0.4 sa_EM 0.05 width 1 end parameters begin compartments EC 3 vol_EC PM 2 sa_PM*width EC CP 3 vol_CP PM EM 2 sa_EM*width CP EN 3 vol_EN EM end compartments
Where a compartment definition syntax is composed of name, dimensions, size and parent compartment. For the purposes of this tutorial it is important to remember that MCell and CellBlender use units of \mu m^3
Hierarchichal View of the Compartments in the Example Model
Once we have defined our compartments, it is possible to define any pattern
that refers to them using one of the two following syntax:
@compartment:species.
For example, let us define the following additional rules to the example model we have been working on:
begin parameters k_r_endo 1 k_recycle 0.1 kp_LL 0.1 km_LL 1 end parameters begin reaction rules ... L(d) + L(d) <-> L(d!1).L(d!1) kp_LL,km_LL @PM:R().R() -> @EM:R().R() k_r_endo @EM:R() -> @PM:R() k_recycle reaction rules
The first rule is a standard rule that specifies that two ligands can bind into a single complex. In cBNGL, rules that do not contain any compartment information imply that they can take place in any compartment where two ligands can be found (following BioNetGen’s don’t show-don’t care principle).
The second rule specifies that whenever it finds two receptors in the plasma
membrane, bound together in some non-specified way, they will be transported to
the endosomal membrane at a k_r_endo rate. Likewise, we specify a
reciclying reaction whereas any kind of receptors can go back to the plasma
membrane. Observables are defined in a similar fashion.
begin observables Species R_Dimers_PM @PM:R.R Species R_Dimers_EM @EM:R.R end observables
In this block we are introducing the Species keyword. This means that, as
we described earlier, instead of counting the number of molecules where a given
pattern appears, it will count the total number of complexes (or species) and
report that to the user.
Let us simulate the system as we have defined it so far: Ligand-receptor binding, ligand-ligand dimerization, receptor phosporylation and receptor transport. Change the simulation so that it runs ten seconds.
Time series of the example model after including compartment information
As we can appreciate, all elements of this simple cascade are being activated and equilibrated (the most downstream element, endosomal receptor, reaches equilibrum at about 8 seconds).
- A pre-generated spatial geometry and reaction networks encoded in the SBML format. It is possible to generate SBML files from BNGL models using the writeSBML() action at the end of your model. If you add this command after a simulate action the resulting SBML model's initial conditions will correspond to the state of your model after the simulation finishes.
- CellBlender
- Plotting software
- full_lr.bngl (Linux) lr_full.blend (MacOSX, Win) (within bngFiles.zip)
SBML is a modeling exchange standard used to encode modeling information in an XML-compliant format. Multiple modeling applications support exporting information as an SBML model (including BioNetGen). Moreover, recent extensions to the SBML standard like SBML-spatial support the definition of 3D geometries. We will make use of these capabilities to define a model in CellBlender
External model import menu
To export an SBML from a BioNetGen file, it is sufficient to add the
writeSBML() action to the bottom of your BioNetGen file.
To import an SBML and SBML-spatial files open Blender and click File, Import, Import External Model. For the purposes of this tutorial you can also use one of the two geometries provided in bngFiles.zip. (icogeometry.mdl and cubicgeometry.mdl) Depending on the size of your model this may take a minute to load. Once loaded you will see imported geometries. The "CellBlender-" sections in the Scene tab on the right Blender panel have also been automatically filled in. As of the time of this tutorial this capability is only supported in Linux.
Before you can simulate the biochemical system you must check that your geometries have consistent normals, are manifold and watertight. Refer to this tutorial for more information on how to make sure your geometry is MCell compliant
Although it is now possible to run your biochemical spatial simulation, you probably want to first set up partitions. These can increase simulation speed by greater than a thousand times! Set your partitions to be 0.1 wide.
The last step before running our simulation is to define the number of iterations our simulation will run for. In the Model Initialization panel increase the number of iterations to 10000.
Save your project to a new blend project file. The model’s reaction output
data will be saved to a folder called projectName_files inside the
directory where the *.blend file is located.
Now we are ready to simulate our file. Click on the run simulation panel and start the simulation.
In order to visualize your data from inside Blender it is sufficient to check the ’Visualize all’ option in the corresponding CellBlender panel. However, given the number of molecules in the example system it is unfeasible to visually store more than a few thousands of iterations. Of more interest is to analyze the behavior of the time series produced by the system; we will explain how to do this in the following section.
It is possible to use any plotting software to visualize the time series generated by MCell/CellBlender, including CellBlender own plugins if you so desire. We have provided a small Python script that satisfies our needs. Running it in the directory where your react_data is located will generate an image for each file found in the directory.
Goal: Be able to generate synthetic geometries and SBML-spatial files from using models learned from fluorescence microscopy images.
MATLAB2013a, CellOrganizer 2.1
Download MATLAB2013a or later and activate it using the installation wizard.
Download and unzip the latest stable release of CellOrganizer 2.1.
Open MATLAB and add navigate to the folder containing the cellorganizer source code by clicking the magnifying glass on the right of the screen. Once in this folder the cellorganizer directory will be visible on the left of the screen in the "CurrentFolder" block. Type \tt{setup} in the Command Window.
Download and extract the SampleData zip file. This folder contains abreviated datasets for the puroposes of demonstration. Place this folder in the same directory as the "cellorganizer" parent folder. Add the "SampleData" to your path by right clicking it and selecting "Add to Path"–>"Selected Folders and Subfolders".
Related CellOrganizer demos: Train, demo2D01, demo3D11, demo3D12, demo3D18, demo3D20
Objectives: To teach the basics of how to train generative models from fluorescence microscopy data.
Procedure: For this tutorial we are going to train a model using 2D images of HeLa cells. These cells have been tagged with a DNA marker, a cytoplasmic marker, and a LAMP2 marker which localizes in lysosomes as seen in :num:`Fig. #sampleimg` below.
LAMP2 tagged HeLa cell
Open the template Train.m by typing "edit Train" in the Command Window block. A file will open and you will see the following
function Train( dnapath,cellpath,protpath,croppath,resolution,filename,dimensionality)
There are clearly several inputs to fill out before you can train a model so let’s go through them in turn.
- The first thing you need to do is tell CellOrganizer where your images are. To
- do this specify the paths as strings like the examples below:
>> dnapath = ['./myDNAimgs/cell*.tif']; >> cellpath = ['./myCELLimgs/cell*.tif']; >> protpath = ['.myPROTimgs/cell*.tif'];
These paths point to the DNA, cytoplasmic, and lysosomal images respectively for a set of cells. The ’*’ wildcard allows you to specify paths with partial patterns.
Warning
If you do not have consistent file names between folders/patterns you may want to list out each file. This will prevent files from becoming mismatched. To do this list the paths in a cell array.
>>dnapath = {'./myDNAimgs/cell1.tif',...
'./myDNAimgs/cell2.tif',...
'./myDNAimgs/cell3.tif'};
It is often useful to specify a cropped region of an image on which to train to avoid noise, artifacts or other cells in different parts of an image that you are not interested in such as the partial cells visible in the bottom of :num:`Fig. #sampleimg`. To do this you will need to specify the croppath that points to a binary mask of your region of interest as well.
>>croppath = ['./crop/cell*.tif']; (If you do not have a mask image, simply enter croppath = [])
You must next specify the resolution at which these images were taken. This is a key part of the model when we want to combine models learned at different resolutions. For the dataset provided in the tutorial, the
resolution = [0.05,0.05]} microns/pixel.
The filename input is a string pointing to where you wish to save the model for example:
>>filename = './tutorialModel';
For the purposes of this tutorial you will be training a 2D generative model on
three cells due to time constraints, so set dimensionality = '2D'. In the
future, to perform 3D training simply set {dimensionality = '3D'} and
change the filepaths to a 3D dataset(see demo3D11).
You are now ready to train a generative model. Use the command
>>Train( dnapath,cellpath,protpath,croppath,resolution,filename,dimensionality)
Additonally a diffeomorphic model may be trained with the inclusion of a boolean variable
>>Train( dnapath,cellpath,protpath,croppath,resolution,filename,dimensionality, isdiffeomorphic)
If you are using the data provided the training should take about an hour.
Note
Larger datasets and 3D data may take hours to train a model so running on a cluster is recommended.**
Once you have started training your model, you may want to explore some of the other training options available in CellOrganizer. Below are a few of the more commonly modified options available. If you modify the template file for your specific uses it is recommended that you save a new version of the file for your specific case.
The following options are editable via the param structure in CellOrganizer. downsampling - This is a 1xD vector where D is the number of dimensions (2 or 3). These numbers are the number of times smaller your images will be for training. This comes at the cost of accuracy for small objects, but can vastly speed up the training process. The default values are [1,1] and [5,5,1] for 2 and 3D training respectively.cytonuclearflag - This is a string flag that may be assigned to ’cyto’,’nuc’,or ’all’ and determines where the vesicular objects are allowed to exist. This parameter is defaulted to ’cyto’ and is crucial to change when working with nuclear proteins.train.flag - This parameter is again a string assigned to ’nuclear’,’framework’, or ’all’. If set to ’nuclear’, only the nuclear model will be trained. Framework will train the nuclear and cell model and ’all’ will train the framework plus a vesicular model.
CellOrganizer dependency structure
Related CellOrganizer demos: Synthesis, demo2D00, demo3D01, demo3D05, demo3D15, demo3DMultiresSynth, demo3DObjectAvoidance, demo3DDiffeoSynth_gmm, demo3DDiffeoSynth_grid, demo3DDiffeoSynth_grid_pick, demo3DDiffeoSynth_uniform.
Objectives To teach the basics of to synthesize images from generative instances from models learned from imaging data.
ProcedureUsing the model created in the previous section you will now
synthesize in silico cells. To do this you will use the Synthesis template.
Enter edit Synthesis into the Command Window to open the template. The file
will read as follows
function Synthesis(modelPath,savePath,numSynthImgs)
If you are using the model created in the previous section set the modelpath input to:
>>modelpath = {'./tutorialModel.mat'};
Next, set the location you would like to save the results. Again, this is just
a string for example savepath = '.'; will save the resulting instances in
the current directory.
Lastly you must set the number of images you would like to sample from the
model. To start you may want to set numSynthImgs = 1;.
Now you are ready to generate an instance from your generative model. Type:
>>Synthesis(modelPath,savePath,numSynthImgs);
This should take a minute or two and generate an image from your model that is saved in the location you specified. This should produce a set of 3D tif images that resemble the slice in figure :num:`Fig. #slice` below.
Single slice from generated LAMP2 pattern
To create an instance containing multiple protein patterns we will simply assign multiple models to the modelpath. For example:
modelpath = {'./models/model1.mat','../models/model2.mat'};
Note
The synthesized models will use the cell and nuclear models contained in the first model in the model list.**
By default, instances from diffeomorphic models are sampled according to the
approximate probability density of the cells and can be synthesized with the
above Synthesis(...) command. In addition to the implicit non-parametric
density sampling, CellOrganizer offers demos to illustrate other sampling
methods, each of which demonstrate different methods of representing the shape
space.
>>demo3DDiffeoSynth_gmm
Marginalizes out the all but the first two dimensions of the shape space and learns a gaussian mixture model over the positions of training images and samples points from this parametric representation.
>>demo3DDiffeoSynth_grid
Synthesizes images from a grid over the first three dimensions of the shape space.
>>demo3DDiffeoSynth_grid_pick
Displays to the user a two-dimensional representation of the shape space and allows the user to select a points to be synthesized.
>>demo3DDiffeoSynth_uniform
Selects a point uniformly at random from the complete shape space.
Once you have started synthesizing your model, you may want to explore some of the other synthesis options available in CellOrganizer. Below are a few of the more commonly modified options available. If you modify the template file for your specific uses it is recommended that you save a new version of the file for your specific case.
The following options are editable via the param structure in CellOrganizer. synthesis - This parameter is string assigned to ’nuclear’,’framework’, or ’all’. If set to ’nuclear’, only the nuclear model will be trained. Framework will train the nuclear and cell model and ’all’ will train the framework plus a vesicular model.output - This is a structure that contains several flags for different types of outputs including ’SBML’,’tifimages’,and ’indexedimage’.
Related CellOrganizer demos: demo3DSBML, demo3DPrimitives, demo3D13,
Objectives To learn how to use an SBML file to automatically generate the necessary geometries for a realistic cellular simulation.
ProcedureTo create a SBML-spatial biochemical system for the SBML file
generated during part two, you will use demo3DSBMLTutorial. Open this file by
typing edit demo3DSBMLTutorial in the Command Window and you will see:
Diffeomorphic shape space with implicit distribution modeled as a gaussian mixture model with 2 components. Color indicates relative probability density, black circles indicate training images in the shape space. In the non-parametric density representation, each simplex (triangle) contains the probability mass.
demo3DSBMLTutorial(SBMLpath,coreSTD,renderSTD).
This demo will read in an SBML file like the one you created in the previous section, attempt to find models that correspond to the compartments in the SBML file, and synthesize images from the appropriate models. This is a modified version of the demo3DSBML file included in the CellOrganizer release. It has been modified to allow you to run the model in a reasonable time.
To use this demo, set the
SBMLpath = './SampleData/motivating_example/Motivating_example_cBNGL2_sbml.xml';
The other two parameters determine how much overlap you will allow the objects
to have and how big they will be respectively. If the coreSTD is set equal to
the renderSTD the objects will not overlap at all, but synthesis will be
dramatically slowed. For demonstration purposes set the coreSTD = 0.5; and
the {renderSTD = 1; to allow the simulation to run quickly. To ensure no
overlap at all, set coreSTD = renderSTD;.
To automatically select the models and generate an SBML-spatial+SBML instance. Enter:
demo3DSBMLTutorial(SBMLpath,coreSTD,renderSTD)
Note
The standard demo3DSBML synthesis will take much longer since it synthesizes high resolution 3D images with no object overlapping allowed.
Sometimes it is very useful to create an SBML-spatial file of just a
pre-synthesized geometry. For example, you may want to use the synthesized
geometries in the provided SampleData folder to create framework meshes. To
do this you will use demo3DImg2SBML. This demo takes in a cell array of
images you wish to create SBML-spatial meshes for. It also takes a resolution
at which the images were synthesized at.
Objectives Now that you’ve become familier with all the tools used to create your simulations run during Part 2, we will analyze the results of these data.
To plot the results of your MCell simulations we will use the {plotMCell}
function. This function takes only one argument {datapath} which is a
string pointing to the ’react_data/seed*/’ folder created by MCell. Where
"seed*" is whatever seed you wish to analyze(e.g. seed_00001). This function will find all the reaction data you have written out from MCell for the
selected seed, create and save a count vs time plot for each species. These
plots will look like :num:`Fig. #lboundem` below.
MCell output of the number of phosphorilated transcription factors vs time (s)
To compare multiple MCell simulations for a specific species, assuming that the simulations have the same outputs (as in the case of this tutorial) use the plotMCellCompare function. This function takes three arguments. The first datapath is similar to the plotMCell function but is a cell array of strings pointing to the paths you wish to compare. For example here we may wish to specify \{datapath = 'cell1/react\_data/seed0001\/','cell2/react\_data/seed0001\/'\}. The second argument here is the species you wish to compare as a number in the file list. For example, in this tutorial Phos_TF would be file number 19. The resulting figure is a plot of the requested species for all the cells requested as seen in :num: Fig. #phosTFcomp below.
Comparison of the phosphorilation of transcription factors for four cells from throughout a shape-space of HeLa cells.
To visualize these results with respect to their location in shape-space you can use the function showCurrentCells(positionpath,imgpath,currcells). Here positionpath specifies the path to the 'selectedPositions.mat' file generated by demo3DDiffeoSynth_grid_pick. Similarly, imgpath is the path to the output images from this demo using a wildcard string for example: './demo3DDiffeoSynth_grid_pick*'. Lastly, currcells specifies the cell numbers you currently have data for or that you plotted using plotMCellCompare. Here you should take care to specify these cells in the same order as you did in plotMCellCompare to ensure the colors across the graphs line up. The resulting figure is a shapespace of the original cells in grey, the specified cells in the color corresponding to their result from plotMCellCompare and the remaining synthetic cells for which you are currently not comparing data in yellow as seen in :num: Fig. #shapeSpaceComp below.
Cells corresponding to the simulated cells in :num: Fig #phosTFcomp allow easy comparison for trends in cellular response across the shape-space.
| [1] | Harris, Leonard A., Justin S. Hogg, and James R. Faeder. "Compartmental rule-based modeling of biochemical systems." Winter Simulation Conference. Winter Simulation Conference, 2009. |
| [2] | T. Peng, Wei Wang, G. K. Rohde1, R. F. Murphy (2009) Instance-Based Generative Biological Shape Modeling. Proceedings of the 2009 IEEE International Symposium on Biomedical Imaging (ISBI 2009), pp. 690-693. |
| [3] | T. Peng and R.F. Murphy (2011) Image-derived, Three-dimensional Generative Models of Cellular Organization. Cytometry Part A 79A:383-391. |