Spatial Stochastic Simulators

1
Spatial Stochastic Simulators

• Kim Avrama Blackwell
• George Mason University
• Krasnow Institute of Advanced Studies

2
Diverse 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
3
Spatial 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
4
MCell
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
5
MCell 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)

 
6
Mcell 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/

7
MCell Model
NMDAR
Pre-synaptic Terminal
Spine Head
Spine Neck
Dendrite
VDCC
Calcium binding Proteins (cytosol)
Pumps (Membrane)
8
UnpairedStimuli

• Calcium differs due to channel distribution

Keller et al. PLoS One 2008
9
Paired Stimuli
EPSP-AP
AP-EPSP

• Calcium depends on timing of AP versus glutamate
  release

Keller et al. PLoS One 2008
10
CDS

• 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

11
CDS 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

12
Stochastic (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

13
Extensions 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

14
Hybrid 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

15
STEPS

• 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

16
STEPS-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
17
Single 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
18
Model 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

19
NeuroRD

• 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

20
NeuroRD 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

21
NeuroRD
22
Calculate j reacting or k moving using Poisson
distribution
23
(No Transcript)
24
NeuroRD - 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
25
NeuroRD

• 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
26
NeuroRD DevelopmentBiochemical Oscillator
Srivastava et al., J Chem Phys
27
Spatial 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)

28
Spatial 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)
29
Spatial Biochemical Oscillator
mRNA
DNA
30
NeuroRD

• 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
31
NeuroRD 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

32
NeuroRD 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

33
NeuroRD 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

34
NeuroRD 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

35
NeuroRD 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

36
NeuroRD 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)

37
NeuroRD 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

38
NeuroRD 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

39
NeuroRD 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