1 / 53

NeuroRD Version 3 July 2016

NeuroRD Version 3 July 2016. Kim “Avrama” Blackwell Zbigniew “Zbyszek” Jedrzejewski-Szmek George Mason University Krasnow Institute of Advanced Studies https://github.com/neurord/stochdiff. Outline. Brief over of algorithm Information required to build signaling pathway models

wildman
Download Presentation

NeuroRD Version 3 July 2016

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. NeuroRD Version 3 July 2016 Kim “Avrama” Blackwell Zbigniew “Zbyszek” Jedrzejewski-Szmek George Mason University Krasnow Institute of Advanced Studies https://github.com/neurord/stochdiff

  2. Outline • Brief over of algorithm • Information required to build signaling pathway models • Experimental design and simulation experiments • Creating model files • Running simulations, and analyzing and visualizing results

  3. NeuroRD • Computationally efficient algorithm to solve the reaction-diffusion master equation (RDME), which is a spatial extension to the chemical master equation • RDME is used to simulate cell signaling pathways, which are biochemical reaction networks in inhomogeneous space • Other RDME algorithms implement an exact stochastic simulation algorithm (SSA) • Steps • MesoRD • StochSS/URDME • Lattice microbes

  4. 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 • Designed for rapid, approximate simulation of large scale models • Numerous reactions and compartments • Large numbers of molecules

  5. NeuroRD • Subdivide dendrites and spines into sub-volumes • Pre-calculate the probability that one molecule leaves the compartment or reacts • Look-up tables store the (binomial) probability that j out of N molecules leave a compartment or react • Original NeuroRD: At each time step, for each molecule, for each voxel: • choose a random number to determine the number, j, molecules out of N leaving or reacting (j is often 0)

  6. NeuroRD v3 – Mesoscopic Simulator • Solves the issue that some molecule quantities are too small to leap • Waste random numbers that result in 0 or 1 firing • Hybrid of tau leap and exact stochastic simulation algorithm without explicit partitioning between single (exact stochastic simulation) and multiple (leaps) events • Continuum of propensities • Propensities change over time

  7. Life cycle of single event channel • For both reaction and diffusion event types • Single events if population or reaction rate (propensity) is small • Multiple events if population or reaction rate is large 2 A ⇋ B 2 voxels:

  8. NeuroRD v3 Simulator • More efficiently deal with large, stiff systems using adaptive, asynchronous t-leaping • Order reaction (including diffusion) channels by time of next event in priority queue; either schedule a leap or single event Highest priority 2 A ⇋ B 2 voxels: Lowest priority Jedrzejewski-Szmek & Blackwell J ChemPhys 2016

  9. NeuroRDv3 Accuracy Comparable to Exact SSA • 1D domain of length 40 • 100 voxels • N=1000 A • 1000 B in lower left corner Jedrzejewski-Szmek & Blackwell J ChemPhys 2016

  10. NeuroRDv3 is ~3x faster than fixed timestept-leaping • Model of dopamine dependent production of cAMP, with feedback via PKA enhancement of PDE Cumulative number of events • Scales linearly with model size Jedrzejewski-Szmek & Blackwell J ChemPhys 2016

  11. Information required to build signaling pathway models • Molecule Species • Reactions and reaction rates • Morphology / shape of space containing the molecules • Discretization of the space • Molecule quantities • Spatial localization and distribution • Diffusion constants • Transient stimulation • Output

  12. Information required to build signaling pathway models • What molecules are important for the phenomenon I’m try to understand? • How does calcium control synaptic plasticity? • Calcium, calcium buffers, calcium pumps, molecules activated by calcium • How do G proteins coupled to nicotinic acetylcholine receptors control growth cones? • Gq subtype of GTP binding protein, Phospholipase C, diacylglycerol

  13. Reaction information • What reactions do these molecules participate in? • Direct reactions: Ca + calbindin <-> CaCalbindin • Indirect reactions (best to avoid these): Camkinase leads to activation of RhoA and Rac 1 (but no enzyme assay has demonstrated direct phosphorylation)

  14. Reaction information • Where to find reaction kinetic values? • Journal of Biological Chemistry, Biochemical Journal, or other modeling papers • Don’t necessarily use modeling papers rate constants, but good source of references • Steady state values and affinities are easier to measure (and find). • Km for enzyme = (Kcat+Kb)/kf • Dissociation constant for bimolecular reactions = Kb/kf • Live cell imaging techniques provide constraint for speed of a small set of reactions • cAMP dynamics in response to dopamine application

  15. Shape of domain (Morphology File) • What is the shape of the cell (or part of cell) I’m trying to model • Dendrite – cylinder of length and radius • Cell body – sphere of radius • Dendrite with spines • Mushroom spines – narrow cylinder+fast cylinder • Growth cone • Large cylinder with several small “finger” cylinders • Use as small and simple a morphology as possible to speed simulations

  16. Quantity (Initial Conditions) • What is the quantity of each molecule in the simulation? • How much ATP in a cell? • How much phospholipase C in the growth cone • What is the location of the molecules • Is the molecule membrane bound or associated or does it diffuse freely? • What is diffusion constant? Can be estimated from molecule weight.

  17. Experimental design • How are you perturbing the tissue? • Current injection • Application of receptor ligand • What are your control simulation experiments? • Compare receptor ligand to no drug control • How does temporal or spatial pattern influence results?

  18. What molecules to evaluate • Possible to evaluate every single molecule in a simulation • Alternatively, only evaluate a subset of molecules • What is the sampling frequency? • High frequency produces huge output files • How quickly do things change?

  19. Model Design Principals • Not all parameters are constrained by literature • Unknown reaction kinetics • Unknown molecule quantities • Model tuning: Adjust unknowns, within reason, to reproduce several sets of data • Model validation: compare to new set of data (not used to tune model) • Make model prediction, follow-up with new experiment

  20. Model Robustness • Parameter sensitivity • How sensitive (or robust) are model results to variation in parameters • Results not same the same as constraining values • E.g. Lowering phosphodiesterase quantity will increase cAMP basal level (constraint). Compensate by increase adenylyl cyclase level. Yields: two good models • A: high PDE, high AC, basal cAMP = 50 nM • B: low PDE, low AC, basal cAMP = 50 nM • Use models A and B for stimulation experiments. • Result: Is peak cAMP the same in both models?

  21. NeuroRD • Model specification allows good experimental design, with separate files for • Reactions (and species specification) • Spatial morphology • Initial conditions • Stimulation • Output specification • Top level file which specifies above files, plus simulation time, spatial grid, random seed, numerical method Tissue Experiment Simulation control

  22. NeuroRD – simulation experiments • Reaction file experiments • What happens if enzyme X is inhibited? • Set kcat = 0 for enzyme X • Compare simulations with and without kcat=0 • Initial Condition file experiments • What happens if knock down enzyme Y? • Set nanomolarity = 0 for enzyme • Compare to simulations with nanomolarity > 0 • Alternative method to inhibit enzyme: • Create reaction of enzyme + inhibitor <-> EI • For inhibition expers, inject inhibitor (change stim file)

  23. NeuroRD – simulation experiments • Morphology file experiments • How is CamKII activation affected by spine morphology? • Create several morphology files with different diameter of spine neck • Compare simulations with different diameters • Stimulation file experiments • How does spacing of spines modify PKA? • Create dendrite with multiple spines • Simulate either two adjacent or two separated spines

  24. NeuroRD – Morphology File • Specify start and end of each segment • Specification includes id, region type, location (x,y,z), radius, and optional label • <Segment id="seg1" region="dendrite"> • <start x="1.0" y="1.0" z="0.0" r="5.0" /> • <end x="1.0" y="11.0" z="0.0" r="5.0" label="pointA"/> • </Segment> • Additional segments start on a previous segment • Id is unique, but multiple segments can be part of same region • Endpoints can be labeled

  25. NeuroRD – Morphology File • Multiple regions • <Segment id="seg1" region="nuc"> • <start x="0.0" y="0.0" z="0.0" r="1" /> • <end x="2.0" y="0.0" z="0.0" r="1" label="pointA"/> • </Segment> • <Segment id="seg2" region="cyt"> • <start on="seg1" at="end"/> • <end x="4.0" y="0.0" z="0.0" r="1" /> • </Segment> • <Segment id="seg3" region="cyt"> • <start on="seg2" at="end"/> • <end x="6.0" y="0.0" z="0.0" r="1" /> • </Segment>

  26. Morphology- Branching • Additional segments connected using “start on=“segment name”, and at=“begin or end” • <Segment id="seg1" region="dendrite"> • <start x="1.0" y="1.0" z="0.0" r="0.5"/> • <end x="3.0" y="1.0" z="0.0" r="0.5"/> • </Segment> • <Segment id="seg2" region="branch1"> • <start on="seg1" at="end" /> • <end x="5.0" y="3.0" z="0.0" r="0.3"/> • </Segment> • <Segment id="seg3" region="branch2"> • <start on="seg1" at="end" /> • <end x="4.0" y="-2.0" z="0.0" r="0.3"/> • </Segment>

  27. Morphology - spines • Randomly distribute spines on dendrite • Shape of spine prototype (multiple allowed) • <SpineType id="mushroom"> • <Section width="0.2" at="0.0" regionClass="neck"/> • <Section width="0.2" at="0.2" regionClass="neck"/> • <Section width="0.2" at="0.3" regionClass="head"/> • <Section width="0.6" at="0.4" regionClass="head"/> • <Section width="0.6" at="0.5" regionClass="PSD"/> • <Section width="0.2" at="0.6" label="pointA"/> • </SpineType> • Distribution • <SpineAllocation id="saA" spineType="mushroom" region="dendrite" lengthDensity="0.6"/>

  28. Morphology • Single spine • Multiple spines of two types • Randomly distributed with specified density

  29. NeuroRD – Reaction File • Define each species • Include diffusion constant, which can be 0 • <Specie name="mGluR" id="mGluR" kdiff="0" kdiffunit = "mu2/s"/> • <Specie name=“cAMP" id=“cAMP" kdiff=“300" kdiffunit = "mu2/s"/> • Specify Reactions • First order – single reactant and product • Second order – two reactants or two products

  30. NeuroRD – Reaction File • Include forward and backward rate constants • If not specified, name and id are the same • <Reaction name = "glu+mGluR--glu-mGluRreac" id="glu+mGluR--glumGluR_id"> • <Reactant specieID="glu" /> • <Reactant specieID="mGluR" /> • <Product specieID="glu-mGluR" /> • <forwardRate> 5e-03 </forwardRate> • <reverseRate> 50e-03 </reverseRate> • </Reaction>

  31. Enzyme reactions • Enzyme reaction is cascade of two reactions, enzyme regenerated • <Reaction name = "DaD1R+Gs--DaD1RGs" id="DaD1R+Gs--DaD1RGs"> • <Reactant specieID="DaD1R" /> • <Reactant specieID="Gsabg" /> • <Product specieID="DaD1RGs" /> • <forwardRate> 0.03e-03 </forwardRate> • <reverseRate> 0.4e-03 </reverseRate> • </Reaction> • 3 products of enzyme reaction: • <Reaction name = "DaD1RGs-DaD1R+GsaGTP+Gbg reac" > • <Reactant specieID="DaD1RGs" /> • <Product specieID="DaD1R" /> • <Product specieID="GsaGTP" /> • <Product specieID="Gbg" /> • <forwardRate> 0.25e-03 </forwardRate> • <reverseRate> 0.0 </reverseRate> • </Reaction> Enzyme Regenerated

  32. NeuroRD – Higher order reactions • Two molecules of same species interact • power=“3” • Stoichiometry of 3, and kinetics from Kf*[CKCamCa4]^3 • <Reaction name = "CKCam bind" id="CKCam_bind"> • <Reactant specieID="CKCamCa4" power="3" /> • <Product specieID="CKCamCa4" power="2" /> • <Product specieID="CKpCamCa4" /> • <forwardRate> 2.0e-12 </forwardRate> • <reverseRate> 0e-3 </reverseRate> • </Reaction>

  33. NeuroRD – higher order reactions • Two molecules bind together, but with 1st order kinetics • n=“2” • Stoichiometry of 2, but kinetics determined from Kf*[cAMP] • <Reaction name = "PKA_bind" id="PKA_bind"> • <Reactant specieID="PKA"/> • <Reactant specieID="cAMP" n="2"/> • <Product specieID="PKAcAMP2"/> • <forwardRate>0.26e-06</forwardRate> • <reverseRate>0.060e-03</reverseRate> • </Reaction>

  34. NeuroRD – Initial Condition File • Initial condition files is dependent on Reaction file, and possibly morphology file (for region specific conc) • Four types of initial conditions • General concentration of molecule in entire morphology: • <NanoMolarityspecieID="mGluR" value="5e3" /> • Region specific concentration • <ConcentrationSet region="PSD" > • followed by <NanoMolarityspecieID=“IP3" value=“30" /> • Overrides general concentration

  35. NeuroRD – Initial Condition File • Surface Density of membrane molecules • Overrides concentration specifications • <PicoSDspecieID="PLC" value="2.5" /> • Region specific Surface Density

  36. NeuroRD – Stimulation File • Stimulation used to inject molecules • Specify molecule species and injection site • Site defined in morphology file • <InjectionStimspecieID="Ca" injectionSite="pointA"> • To inject at several spines, called sa1: • <InjectionStimspecieID="a" injectionSite="sa1[3,4,5].pointA"> • To inject at all spines, called sa1: • <InjectionStimspecieID="a" injectionSite="sa1[:].pointA"> • To inject at submembrane of region “dend”: • <InjectionStimspecieID="a" injectionSite=“dend:submembrane">

  37. NeuroRD – Stimulation File • Repetitive trains can be created • Specify onset time, duration, rate (amplitude) • Optional: period and end used for train • Optional: InterTrain Interval, numTrains to repeat train onset end duration rate InterTrainInterval period

  38. Stimulation File • <InjectionStim specieID="glu" injectionSite="sa1[:].pointA"> • <onset> 1000 </onset> • <duration> 30 </duration> • <rate> 40 </rate> • <period> 100 </period> • <end> 2000 </end> • <interTrainInterval> 4000 </interTrainInterval> • <numTrains> 2 </numTrains> • </InjectionStim>

  39. Arbitrary input patterns • Injection specified with filename containing time and rate: • <InjectionStimspecieID=“calcium" injectionSite="pointA"> • <rates> • <xi:includehref="rates.txt" parse="text" /> • </rates> • </InjectionStim> • Rates can be outputs from electrical model, e.g. genesis, neuron or moose • Rates can be arbitrary patterns, format: • Time1 value • Time2 value • …

  40. NeuroRD – Output Specification • Specify outputset name (filename), dt for output, species • Optionally specify region, e.g. dendrite • <OutputSet filename = “denddt1" region="dendrite" dt="1.0"> • <OutputSpecie name="glu" /> • <OutputSpecie name="IP3" /> • </OutputSet> • Multiple outputSets can be specified • Sample slowly changing molecules less frequently • <outputSet filename=“dt10” dt=“10.0”> • Sample glutamate receptors from PSD only • <Outputset filename=“psd” region=“PSD” dt=1.0>

  41. NeuroRD version 3 – Model file • Specify all the other files: • <?xml version="1.0" encoding="UTF-8" standalone="yes"?> • <SDRunxmlns:xi="http://www.w3.org/2001/XInclude" xmlns="http://stochdiff.textensor.org"> • <xi:includehref="Purkreactions " /> • <xi:includehref="Purkmorph " /> • <xi:includehref="Purkstim " /> • <xi:includehref="Purkic " /> • <xi:includehref="Purkio " />

  42. NeuroRD - Model file • Algorithm for simulations: • New (faster, more accurate) algorithm <calculation>GRID_ADAPTIVE</calculation> • accuracy control parameter: <tolerance> 0.01 </tolerance> • Original, fixed time step tau leap: • <calculation>GRID_STEPPED_STOCHASTIC <calculation> • time step : <fixedStepDt> 0.005 </fixedStepDt> • Deterministic: • <calculation>GRID_STEPPED_CONTINUOUS<calculation> • random seed – change stochastic realization • Simulation seed • Spine seed (for distribution of spines)

  43. NeuroRD - Model file • Indicate total simulation time (in msec) • <runtime> 300000 </runtime> • Interval (infrequent) for writing output of every molecules in every voxel • can be used to re-start simulations, or for initial conditions of other simulations. • <outputInterval> 500.0 </outputInterval> • Geometry can be 2D: 1 layer of voxels of specified depth: • <geometry> 2D </geometry> • <depth2D> 0.6 </depth2D> • Geometry can be 3D: cylindrical mesh to maintain computational efficiency

  44. Effect of 2D depth • Smaller depth has smaller volume and noisier molecule concentrations

  45. 3D morphology • Cylindrical elements allow larger and fewer mesh elements. • Update with 2D vs 3D comparison of branched morphology (or branch with spines)

  46. Discretization (mesh or voxels) • General, or region specific • Spines are specified separately – 1D diffusion • Smaller submembrane allowed • Subdivide submembrane voxels into smaller voxels using <surface layers> • <discretization> <defaultMaxElementSide>0.3</defaultMaxElementSide> • <spineDeltaX> 0.1 </spineDeltaX> • <surfaceLayers>0.1,0.2</surfaceLayers> • <MaxElementSide region=“soma">1.0</MaxElementSide> </discretization>

  47. Mesh element size • Effect of surface layers: subdivide the submembrane voxels into smaller voxels: • Specify coarser mesh size for large parts:

  48. NeuroRD – running simulation • Java -jar neurord-3.1.3-all-deps.jar Model.xml • Current users guide (readme): github.com/neurord/stochdiff • Latest release: github.com/neurord/stochdiff/releases • h5 output file is default • Python programs for Output analysis • github.com/neurord/D1pathways/: • neurordh5_analysis.py, h5utils.py, plot_h5.py • Create output files with region specific time traces • Plot overall concentration versus time • Calculate a basal value, peak, min, peaktime, etc.

  49. NeuroRDViz • Visualization program written in python • Most morphology pictures in this presentation created using it • Movies of change in concentration over time and space

  50. NeuroRDViz • Visualization program written in python • Most morphology pictures in this presentation created using it • Movies of change in concentration over time and space • Demo

More Related