Title: Spatial Stochastic Simulators
1Spatial Stochastic Simulators
- Kim Avrama Blackwell
- George Mason University
- Krasnow Institute of Advanced Studies
2Diverse Numbers of MoleculesSpatially
Inhomogeneous
G protein coupled receptors Diffusion required
for signal interaction
Glutamate receptors 1 ?M ? 60 molecules molecular
interactions occur stochastically
Small number of molecules in spines Large number
of molecules in system
Kotaleski and Blackwell 2010
3Spatial Stochastic Simulators
- Particle based
- Smoldyn, MCell, CDS
- Individual molecules are represented as
point-based particles, which diffuse random
distance and random direction at each time step - If two reacting molecules pass near each other
they may react - Computations increase with number of molecules
Diffusion
Association
Dissociation
sb
Membrane
4MCell
Transparent
- Geometry from volumetric imaging data using
Blender (www.blender.org) - Mesh elements may be reflective, transparent, or
absorptive - Surface or volume diffusion
- Ray tracing determines whether molecules would
have collided during (fixed) time step - Reaction rules depend on order of reaction, and
whether surface or volume molecules involved
(Kerr et al. SIAM J Sci Comput 2008)
Reflective
5MCell Diffusion
- Diffusion distance from probability density
- Radial distance from uniformly distributed random
variable X - Speed computations by storing values of X in
look-up table - Direction uniformly distributed random variable
0,2p)
Â
6Mcell STDP example
- Calmodulin activation versus spike timing
- Do NMDA receptors and VDCC produce different
calmodulin profiles? - Neuron model to determine voltage-dependent open
probability of VDCC and NMDA - MCell model with calmodulin, calbindin, NCX and
PMCA - Model Keller et al. PLoS One 2008, tutorial
http//www.mcell.org/tutorials/
7MCell Model
NMDAR
Pre-synaptic Terminal
Spine Head
Spine Neck
Dendrite
VDCC
Calcium binding Proteins (cytosol)
Pumps (Membrane)
8UnpairedStimuli
- Calcium differs due to channel distribution
Keller et al. PLoS One 2008
9Paired Stimuli
EPSP-AP
AP-EPSP
- Calcium depends on timing of AP versus glutamate
release
Keller et al. PLoS One 2008
10CDS
- Particle based simulator with event driven
algorithm - All possible collisions are detected during short
dt - If collision detected, the exact collision time
is calculated - Earliest collision (or reaction events) are
simulated one-by-one until dt - Particles have volume, thus can simulate crowding
and volume exclusion - http//nba.uth.tmc.edu/cds/content/download.htm
11CDS Example
- Morphology from triangular meshes
- CaMKII diffusion out of spine depends on
morphology (b) and also binding targets and
F-actin - Byrne et al. J Comput Neuro 2011
12Stochastic (non-spatial) Simulators
- Gillespie (Exact Stochastic Simulation Algorithm)
Propensity of reaction aj ? Kf ? Np - Propensity of any reaction, a0 ? aj
- Next reaction occurs with exponential
distribution with mean a0 - Identity of reaction selected randomly, based on
propensity - Computations increase with number of molecules
13Extensions to Gillespie Algorithms
- Spatial Gillespie, e.g. Fange et al. 2010, PNAS
- Morphology is subdivided into small compartments
- Propensity of diffusion calculated from diffusion
coefficient, ad ? D ? Nd - Diffusion considered as another reaction
- Tau leap non-spatial
- Allow multiple reaction events, Kj, to occur for
each reaction at each time step, t, according to
Poisson
14Hybrid Models
- Partition the reaction-diffusion space into two
or more sets of reactions (and diffusion) - Each set is simulated differently
- Diffusion deterministic, reactions stochastic
- Fast reactions - deterministic, slow reactions
stochastic - Critical reactions - exact stochastic,
non-critical reactions tau leap
15STEPS
- Spatial extension of exact stochastic simulation
algorithm - Tetrahedral meshes allows realistic geometries
- Diffusion constant can vary between compartments
- Simulations are specified in python, witih
morphology, reactions and simulations specified
independently (for ideal control of simulation
experiments) - http//steps.sourceforge.net/STEPS/Home.html
16STEPS-Cerebellar LTD
Calcium Buffers
Calcium Pumps
AMPA Receptor
calcium
Protein Phosphatase 5
PKC
Raf
Raf-act
Protein Phosphatase 2A
Arachidonic Acid
Positive Feedback Loop
MEK
MapKinase Phosphatase 1
ERK
cPLA2
Protein Phosphatase 1
Inactivation, dephosphorylation Activation,
phosphorylation
17Single Spine Model
- Average of multiple simulations reveals graded
induction of LTD - Single runs reveals bistability at intermediate
calcium
Time (min)
Antunes et al. J Neurosci 2012
18Model Limitations
- All these model have either small volume (single
spine) or small number of reactions
(calmodulinCaMKII) - Only MCell model uses voltage to determine
calcium influx - Smoldyn
- Particle simulation algorithm incorporated into
Moose (Genesis 3) and VCell - No neuroscience examples yet
19NeuroRD
- Spatial extension to Gillespie tau leap
- Multiple reaction events and diffusion events can
occur during each time step - Morphology is subdivided into small compartments
- Cuboidal meshes and cylindrical meshes possible
20NeuroRD Mesoscopic
- Subdivide dendrites and spines into sub-volumes
- Pre-calculate the probability that one molecule
leaves the compartment or reacts
- Look-up tables store the probability that j out
of N molecules leave a compartment or react - At each time step, for each molecule, choose a
random number to determine the number, j,
molecules out of N leaving or reacting
21NeuroRD
22Calculate j reacting or k moving using Poisson
distribution
23(No Transcript)
24NeuroRD - Validation
- An approximation, to allow large scale
simulations - Agrees with Smoldyn, and deterministic solution
for reaction-diffusion system
Molecule A
Molecule B
Oliveira et al. 2010, PLoS One
25NeuroRD
- NeuroRD is up to 60 times faster than Smoldyn
- Computations increase linearly with number of
compartments, but not molecules
NeuroRD NeuroRD Smoldyn Smoldyn
Simulation initial molecules injected Time (hmmss) Memory (kb) Time (hmmss) Memory (kb)
Diffusion 0 2000 00002.86 1608 00007.04 2344
Reaction 28853 0 00005.97 1764 00803.53 26524
Reaction Diffusion I 662 4000 00004.51 1764 00248.90 22168
Reaction Diffusion II 6619 40000 00007.58 1772 21958.00 23760
Oliveira et al. 2010, PLoS One
26NeuroRD DevelopmentBiochemical Oscillator
Srivastava et al., J Chem Phys
27Spatial Gene Oscillator
- mRNA is inactive in the nucleus, diffuses into
cytosol - A diffuses to nucleus, binds to DNA
- Effect of diffusion constant (2 cytosol
compartment)
28Spatial Biochemical Oscillator
Inactive mRNA in nucleus, activated by binding in
cytosol compartment Vary number of compartments,
and translation compartment
Protein quantity
mRNA production is faster when A binds to
DNA mRNA production and degradation are faster
for A than R Protein synthesis and degradation
are faster for A than R R degrades A (at same
catalytic rate that A spontaneously degrades)
29Spatial Biochemical Oscillator
mRNA
DNA
30NeuroRD
- Model specification allows good experimental
design, with separate files for - Reactions
- Spatial morphology
- Initial conditions
- Stimulation
- Output specification
- Top level file which specifies reactions,
morphology, initial conditions, output specs,
time step and spatial grid, random seed
Tissue
Experiment
Simulation control
31NeuroRD Morphology File
- Specify start and end of each segment
- Specification includes id, region type, location
(x,y,z), radius, and optional label - ltSegment id"seg1" region"dendrite"gt
- ltstart x"1.0" y"1.0" z"0.0" r"0.5" /gt
- ltend x"1.0" y"2.0" z"0.0" r"0.5"
label"pointA"/gt - lt/Segmentgt
- Additional segments start on a previous segment
- Branching is possible see branching.tar
32NeuroRD Reaction File
- Define each species that has either a reaction
pool or conservepool - Include diffusion constant, which can be 0
- ltSpecie name"mGluR" id"mGluR" kdiff"0"
kdiffunit "mu2/s"/gt - Specify Reactions
- First order single reactant and product
- Second order two reactants or two products
33NeuroRD Reaction File
- Include forward and backward rate constants
- ltReaction name "glumGluR--glu-mGluR reac"
id"glumGluR--glumGluR_id"gt - ltReactant specieID"glu" /gt
- ltReactant specieID"mGluR" /gt
- ltProduct specieID"glu-mGluR" /gt
- ltforwardRategt 5e-03 lt/forwardRategt
- ltreverseRategt 50e-03 lt/reverseRategt
- ltQ10gt 0.2 lt/Q10gt
- lt/Reactiongt
34NeuroRD Initial Condition File
- Four types of initial conditions
- General concentration of molecule in entire
morphology, or - Region specific concentration
- Overrides general concentration
- Surface Density of membrane molecules
- Overrides concentration specifications
- Surface Density of Membrane molecules in specific
region - Overrides general surface density
35NeuroRD Initial Condition File
- General concentration of each molecule should be
specified (zero otherwise) - ltNanoMolarity specieID"mGluR" value"5e3" /gt
- Surface density if molecule is membrane bound
- ltPicoSD specieID"PLC" value"2.5" /gt
- Initial conditions for different parts of
morphology - ltConcentrationSet region"PSD" gt
- followed by ltNanoMolarity specieIDIP3"
value30" /gt
36NeuroRD Stimulation File
- Stimulation used to inject molecules
- Temporary fix until software is integrated with
software for simulating neuron electrical
activity and ion channels - Specify molecule and injection site
- ltInjectionStim specieID"Ca" injectionSite"pointA
"gt - Repetitive trains can be created
- Specify onset time, duration, rate (amplitude)
- period and end used for train
- InterTrain Interval to repeat train (e.g. For LTP)
37NeuroRD Output Specification
- Specify dt for output, species and compartment
- ltOutputSet filename "dt1" region"dendrite"
dt"1.0"gt - ltOutputSpecie name"glu" /gt
- ltOutputSpecie name"IP3" /gt
- lt/OutputSetgt
- Multiple outputSets can be specified
- Sample slowly changing molecules less frequently
- Sample glutamate receptors from PSD only
38NeuroRD Model file
- Specify all the other files
- ltreactionSchemeFilegtPurkreactionslt/reactionSchemeF
ilegt - ltmorphologyFilegtPurkmorphlt/morphologyFilegt
- ltstimulationFilegtPurkstimlt/stimulationFilegt
- ltinitialConditionsFilegtPurkiclt/initialConditionsFi
legt - ltoutputSchemeFilegtPurkiolt/outputSchemeFilegt
- Specify some other parameters, such as algorithm
variations and random seed - Indicate total simulation time, time step and
largest compartment size
39NeuroRD running simulation
- Java -jar stochdiff.jar Purkmodel.xml
- Morphology output file
- Purkmodel.out-mesh.txt
- Ascii output file
- name of model file -- .out output set name -
conc.txt - Purkmodel.out-dt1-conc.txt