Integration of rulebased models and compartmental models of neurons
Abstract
Synaptic plasticity depends on the interaction between electrical activity in neurons and the synaptic proteome, the collection of over 1000 proteins in the postsynaptic density (PSD) of synapses. To construct models of synaptic plasticity with realistic numbers of proteins, we aim to combine rulebased models of molecular interactions in the synaptic proteome with compartmental models of the electrical activity of neurons. Rulebased models allow interactions between the combinatorially large number of protein complexes in the postsynaptic proteome to be expressed straightforwardly. Simulations of rulebased models are stochastic and thus can deal with the small copy numbers of proteins and complexes in the PSD. Compartmental models of neurons are expressed as systems of coupled ordinary differential equations and solved deterministically. We present an algorithm which incorporates stochastic rulebased models into deterministic compartmental models and demonstrate an implementation (“KappaNEURON”) of this hybrid system using the SpatialKappa and NEURON simulators.
Keywords:
Hybrid stochasticdeterministic simulations, hybrid spatialnonspatial simulations, multiscale simulation, rulebased models, compartmental models, computational neuroscienceM
Authors’ Instructions
1 Introduction
The experimental phenomena of long term potentiation (LTP) and long term depression (LTD) show that synapses can transduce patterns of electrical activity on a timescale of milliseconds in the neurons they connect into longlasting changes in the expression levels of neurotransmitter receptor proteins. This synaptic plasticity plays a crucial role in the development of a functional nervous system and in encoding semantic memories (e.g. motor patterns) and episodic memories (experiences), converting stimuli lasting for seconds into memories that last a lifetime [16].
There are a number of computational models of how synaptic plasticity arises from patterns of pre and postsynaptic electrical activity, the dynamics of amino3hydroxy5methyl4isoxazolepropionic acid receptors (AMPARs) and NmethylDaspartic acid receptors (NMDARs), calcium influx through these receptors and intracellular signalling in the postsynaptic density (PSD), a dense, proteinrich structure attached to neurotransmitter receptors [1, 22, 14, 28]. The level of detail of the molecular component of these models ranges from deterministic simulations in one compartment [1] through stochastic models with coarse granularity [27] and, at the most detailed, particlebased simulations in which the Brownian motion of individual molecules is modelled [26, 28].
The model with the greatest number of molecular species has 75 variables representing the concentrations of signalling molecules, complexes of signalling molecules and phosphorylation states [1]. This constitutes a small subset of the 1000 proteins identified in the mouse postsynaptic proteome, the collection of proteins in the PSD [6]. Even the subset of the postsynaptic proteome containing proteins associated with membranebound neurotransmitter receptors contains over 100 members [20]. As these proteins are particularly associated with synaptic plasticity, it would be desirable to increase the number of proteins and complexes it is possible to describe in simulations.
As each protein has a number of binding sites, a combinatorially large number of complexes can arise, meaning that a correspondingly large number of reactions are needed to describe the dynamics of all possible complexes. By specifying rules whose elements are fragments of complexes, rulebased languages and simulators [5], such as Kappa [7] or BioNetGen [9], obviate the need to specify reactions for all possible complexes. These rules are simulated using a method similar to Gillespie’s stochastic simulation method for reactions [7]. Kappa has been used to predict the sizes of clusters of proteins in the postsynaptic proteome [23].
Compartmental models of electrical activity in neurons split the neuronal morphology into a number (ranging from 1 to around 1000) of compartments, and specify the dynamics of the membrane potential in each compartment in terms of coupled ordinary differential equations (ODEs) [12, 25]. Quantities beyond the membrane potential can also be modelled, e.g. intracellular calcium concentration and concentrations of a few other molecules such as buffers and pumps. Various packages can generate and solve the equations underlying compartmental models from various model description languages, for example NEURON [3], MOOSE [21] and PSICS [2].
We present an algorithm which integrates rulebased models and compartmental models of neurons. To be sure of understanding a simple, yet interesting, case, we limit ourselves to considering isolated postsynaptic proteomes in a neuron of arbitrary morphology. Although of interest, we do not consider diffusion of molecules within the neuron. We implement the algorithm by incorporating the SpatialKappa rulebased simulator [24] into the NEURON simulator [3]. We validate the combined simulator (“KappaNEURON”) against standalone NEURON, and demonstrate how the system can be used to simulate complex models.
2 Simulation Method
2.1 The System to Be Simulated
An example of the type of system to be simulated is shown in Fig. 1. There is a hippocampal CA1 pyramidal neuron (left) upon which are located a number of synapses. Excitatory synapses are generally located on synaptic spines, small protuberances from the neuron whose narrow necks limit, to an extent, diffusion of ions and molecules between the spine head and the rest of the neuron (middle). The synapse contains a postsynaptic proteome of arbitrary complexity (right). Firing in the presynaptic neuron causes release of neurotransmitter from the presynaptic bouton, which, after diffusing across the synaptic cleft, binds to AMPARs and NMDARs. The AMPARs open and close on a submillisecond timescale, allowing sodium ions to flow into the cell. These ions charge the membrane locally and flow to other parts of the neuron, where they also charge the membrane (middle, Fig. 1). The NMDARs open and close with slower dynamics, and allow calcium ions to flow into the spine. Inside the spine, the calcium ions bind to various proteins such as calmodulin and the resulting calciumcalmodulin complex may then bind to Ca/calmodulindependent protein kinase II (CaMKII), initiating signalling known to be critical for the induction of LTP and LTD.
We will first give the general set of equations that constitute a deterministic description of this type of system, then apply the general equations to a specific example, and next describe how the equations are solved. Finally we will show how we incorporate rulebased models and simulate the hybrid system.
2.2 Deterministic Description of Electrical and Chemical Activity in a Neuron
To describe electrical activity in the neuron, the neuron is split into compartments, each of which should be small enough to be approximately isopotential. Apart from the root compartment, which is located in the soma, each compartment has a parent, and a compartment may have one or more children. The equation for the membrane potential in compartment derives from Kirchoff’s current law:
(1) 
The left hand side is the current per unit membrane area charging or discharging the membrane; is the specific membrane capacitance in compartment . The first term on the right hand side describes current flow into compartment from its neighbours ; is the intracellular resistivity, is the path length between the midpoints of and each of its neighbours , and is the mean diameter of the path. The second term on the right hand side is the total transmembrane current per unit area (referred to as current density) in compartment carried by various species of ion via ion channels () and membrane pumps (), which act to maintain concentration differences between the intracellular and extracellular space. To represent “nonspecific” ion currents whose concentration is not accounted for, there is final term . Here the minus sign reflects the conventions that inward current is negative and the extracellular space is regarded as electrical ground.
The current density carried by species through types of ion channel is:
(2) 
where is the conductance of ion channel type in compartment , which may be a function of time or a state variable , and and is a function describing the – characteristic of current flow of ions of type through channel , which may depend on , the intracellular concentration of in compartment , and , the extracellular concentration of , which is assumed to be constant. A normalised GoldmannHodgkinKatz (GHK) current equation [25] can be used for . For channels through which calcium flows, the typically large ratio between intracellular and extracellular calcium concentrations means this function depends quite strongly on the intracellular calcium concentration , but in channels not permeable to calcium it is usual to use a linear approximation , where is the reversal potential for that channel. By removing the dependence on intracellular concentrations , this approximation allows currents carried by ions other than calcium to lumped together in a nonspecific ion category.
The state variable is the number of type ion channels in compartment which are in an open conformation. It is modelled as the occupancy of a state of a Markov process with membrane potentialdependent transition rates. For small number of channels the Markov process is simulated stochastically, but with large numbers of channels the system is practically deterministic and the master equation of the Markov process is simulated using ODEs.
The dynamics of the intracellular concentrations of ions can be modelled using further ODEs. The rate of change of depends on , the channel transmembrane current density carried by , and consumption and release by intracellular reactions:
(3) 
where is the surface area of the compartment, is the volume of the compartment, is the valency of ion , and is Faraday’s constant. The term describes the net flux of due to intracellular reactions . It arises from treating the intracellular reactions in compartment as a set of kinetic schemes:
(4) 
The flux of arising from this reaction would be:
(5) 
The pump current may be defined in terms of the flux of a reaction , for example:
(6) 
where the prefactor converts from flux to current. Thus the whole electrical and molecular system is defined by a system of ODEs.
2.3 Example Deterministic Description
To help understand the formalism above, we provide an example of a simple onecompartment system; this will also be used as the validation example in Section 4.1. There is a single compartment whose membrane contains passive (leak) channels, calcium channels, and a transmembrane calcium pump, described by the kinetic scheme:
Ca binding:  (7)  
Ca release: 
where Ca represents intracellular calcium, P represents a pump molecule in the membrane, \ceP.Ca is the pump molecule bound by calcium and and are rate coefficients.
We substitute and the flux of the “Ca release” reaction into Equation (6) to obtain the pump current:
(8) 
where we have used the fact that the total concentration of the pump molecule is the sum of the concentrations [P] and [\ceP.Ca] of unbound and bound pump molecule. Since there is only a single compartment, we have dropped subscripts. The calcium channel current flowing into the compartment is determined by Equation (2) with a constant conductance of magnitude , and the nonspecific current is used for the passive channels so that , where is the passive reversal potential.
To construct the ODEs corresponding the kinetic scheme (7) and the expression for the pump current (8), Equation (5) is applied to the scheme to give fluxes, which, along with , and , are substituted in Equations (1) and (3) to give:
(9)  
(10)  
(11) 
The notional volume may describe the volume of a thin submembrane shell rather than the volume of the whole compartment. We assume that is the volume of the whole cylindrical compartment so that , where is the diameter of the compartment.
2.4 Simulation of Deterministic Variables
Simulators of deterministic electrical and chemical activity in neurons, such as NEURON [3], solve the coupled ODEs by gathering the variables , and other state variables into one state vector and solving the ODE system:
(12) 
where is the rate of change of each state variable and is a time dependent forcing input. In principle depends on all variables, though the structure of compartmental models means each element of depends on only a few elements of . These equations can be solved by implicit Euler integration, which, although not providing the secondorder accuracy of more advanced schemes, does give guarantees about numerical stability and is used by default in NEURON [18]. In implicit Euler the derivative is evaluated at , the end of the time step:
(13) 
Taylor expanding the right hand side of (13) in and rearranging gives the equation for one update step:
(14) 
where is the Jacobian matrix at time , which can be computed numerically, or analytically for efficiency. To optimise simulation speed, the time step may vary depending on the rate of change of the variables; but whatever the value of , all variables are updated simultaneously.
2.5 Modifications to Accommodate RuleBased Simulation
To combine fixedstep simulation of continuous variables described by ODEs with discrete variables described by stochastic, rulebased models, we use principles akin to those used in hybrid simulations of systems of chemical reactions [13]. The state variables (elements of ) are partitioned into continuous variables, which are updated at fixed intervals of by an ODE solver, and discrete variables, which are updated asynchronously by the rulebased solver, as outlined in Appendix A.1. Some “bridge” variables are referred to by both solvers. The combined simulation algorithm must ensure that the two solvers are synchronised appropriately and that conversions between continuous and discrete quantities are made.
For simulations combining molecular and electrical activity (e.g. Fig. 1) the membrane potential would be a continuous variable, intracellular molecules such as calmodulin and CaMKII would be stochastic variables, and the intracellular calcium in the spine would be a stochastic bridge variable.
Conversions
In deterministic simulations of biochemical reactions in neurons (e.g. [1]) a molecular species or ion is represented by an intensive quantity – its concentration ; whereas in stochastic simulations it is represented by an extensive quantity – the number of molecules in the volume in which exists. Thus to compare the deterministic (ODE) and stochastic (rulebased) parts of the simulation, intensive and extensive quantities need to be interconverted using Avogadro’s constant :
(15) 
Rate coefficients for reactions based on concentrations must also be converted to ones appropriate for species number for use in the rulebased simulator’s rules. To derive the conversion formula, consider a bimolecular kinetic scheme in which is the rate coefficient, i.e. \ceS + T ¿[k] S.T. Typical units for are \reciprocal\molar\reciprocal\second. The kinetic scheme can be expressed as an ODE . Converting the concentrations according to (15) yields an equivalent ODE whose variables are numbers of molecules: , where is the converted rate coefficient and has units \reciprocal\second. In general for an equation with reactants, the relation between rate coefficients for numbers of molecules and concentrations is . State variables (for example the state of channels) may also be controlled by the rulebased simulator, and here other conversion formulae apply.
Creation and Destruction Rules for Bridging Variables
To simulate the channel currents in the rulebased solver we need to write creation or destruction rules that are equivalent to in Equation (3). These rules are:
(16)  
(17) 
Here can be an expression that references continuous or discrete variables.
Update and Synchronisation
The procedure for updating the time from to (Fig. 2) is:

Pass all relevant continuous variables, e.g. conductances and voltages needed to compute in the rulebased simulator.

Run the rulebased simulator from to .

Compute the net change in the total number of each bridging species (including in any complexes) in compartment over the time step and convert each change back to a current: . For each membrane potential , set the corresponding element of equal to (see second term on right of Equation 1).

Update the continuous variables according to the update step (14).
When running the rulebased model, it will not stop precisely on the boundary of the time step since the update times are generated stochastically. To deal with this problem, the time of the next event in the rulebased component is computed before updating the variables. As soon as the next event time is after the end of the deterministic step, that event is thrown away, as justified in Appendix A.2.
3 Implementation
We have implemented the algorithm described in the previous section by
linking the Javabased SpatialKappa implementation of the Kappa
language [24] to version 7.4 of NEURON, which
allows reactiondiffusion equations to be specified in Python
[18]. Our implementation (“KappaNEURON”) is
available at
http://github.com/davidcsterratt/KappaNEURON.
We have used the
py4j
In order to specify a model the Kappa component is specified in a separate file. Fig. 3 shows an example of a simple calcium pump specified in Kappa. This file is then linked into the NEURON simulation as demonstrated in the Python code in Fig. 4. The mechanisms specified in the Kappa file take over all of NEURON’s handling of molecules in the cytoplasm of chosen sections.
4 Results
4.1 Validation
We validated our implementation by comparing the results of simulating ODE and rulebased versions of the model described in Section 2.3 using standard NEURON and KappaNEURON respectively. Fig. 5A shows the deterministic solution of the system of ODEs (9)–(11) (blue) and a sample rulebased solution using the Kappa rules in Fig. 3 (red) from a single compartment with diameter 1\micro\meter and length 1\micro\meter, giving a volume within the range 0.01–1\micro\meter\cubed typical of spine heads in the vertebrate central nervous system [11].
The calcium conductance is zero apart from during a pulse lasting from 5–10\milli\second when an inward calcium current begins to flow ( is negative). This causes a sharp rise in intracellular calcium concentration [Ca], which, because of the GHK current equation, reduces the calcium current slightly, accounting for the initially larger magnitude of the calcium current. As the calcium concentration increases, it starts binding to the pump molecules, depleting the amount of the free pump molecules [P]. Once the calcium channels close (), the calcium influx stops, and the remaining free calcium is taken up quickly by the pumps. The pumpcalcium complex dissociates at a slower rate, leading to a positive (outwards) calcium current. The stochastic traces (red) are very similar to their deterministic counterparts (blue) apart from some random fluctuations, particularly in the trace of calcium. This agreement, along with a suite of simpler tests included with the source code, indicates that the implementation is correct.
Fig. 5B shows deterministic (blue) and stochastic (red) simulations in a spine with a diameter of 0.2\micro\meter (i.e. 1/25 of the volume of the simulation in Fig. 5A). The shape of the traces differs due to the change in surface area to volume ratio. Due to the smaller volume and hence smaller numbers of ions involved, the fluctuations are relatively bigger.
4.2 Demonstration Simulation
To demonstrate the utility of integrated electrical and rulebased neuronal models, we constructed a model of a subset of the synaptic proteome, with a focus on the signal processing at the early stages of the CaMCaMKII pathway. We encoded in Kappa published models of: the dynamics of NMDARs [27]; binding of calcium with calmodulin and binding of calmodulincalcium complexes to CaMKII [19]; and binding of calcium to a calbindin buffer [8]. We embedded these linked models into a simple model of a synaptic spine, comprising head and neck compartments, connected to a dendrite. As well as the NMDARs, modelled in Kappa, there were AMPARs in the spine head, and backpropagating spikes were modelled by inserting standard HodgkinHuxley ion channels [12] in the dendritic membrane. To emulate spiketiming dependent synaptic plasticity protocols [15], a train of 10 excitatory postsynaptic potentials (EPSPs) were induced in the synapse in the spine head at 20Hz, each of which was followed by an action potential. There were also 50 other synaptic inputs onto the dendrite, though these did not contain the rulebased model.
Fig. 6 shows results from one simulation at short (0.6s) and long (6s) time scales. The first stimulation of the detailed synapse paired with action potential initiation occurs at 0.1s, as can be seen by the voltage trace. The stimulation releases glutamate, which binds to AMPARs (which remain open for a few milliseconds) and NMDARs (which remain open for 100s of milliseconds). Due to the backpropagating action potentials releasing the voltagedependent block of NMDARs, there are peaks in the calcium current () at the same time as the action potentials. The calcium entering through the NMDARs binds to the calbindin and calmodulin (CaM) buffers. The CaMCa complex can bind to CaMKII, the rate depending on which of the C and N lobes of CaM the Ca are bound to. This CaMKIICaMCa complex can then be phosphorylated, leading to a long lasting elevation in its level. This will then to phosphorylate stargazin, which will help to anchor AMPARs in the membrane, thus contributing to LTP.
5 Discussion
We have presented a method for integrating a stochastically simulated rulebased model of proteins in a micronsized region of a neuron into a compartmental model of electrical activity in the whole neuron. The rulebased component allows the biochemical interactions between binding sites on proteins to be specified using a tractable number of rules, with the simulator taking on the work of tracking which complexes are present at any point during the simulation.
Our approach is similar to that of Kiehl et al. [13] who simulated chemical reactions with a hybrid scheme. However their integration scheme synchronised at every discrete event, whereas in ours, synchronisation is driven only by the time step of the continuous simulator. This principle is appropriate for neural systems, in which we can expect many discrete events per time step.
Recently Mattioni and Le Novère [17] have integrated the ECell simulator with NEURON. Our approach is similar to theirs, though with two differences. Firstly we have integrated a rulebased simulator. This has the advantage that the interactions between the combinatorially large numbers of complexes present in the PSD can be specified using a tractable number of rules, though this limits the simulation method to the sequential style of Gillespie’s stochastic simulation algorithm, and does not allow for any of the approximations that increase the algorithm’s efficiency. Secondly Mattioni and Le Novère get the ODEbased solver to handle calcium, whereas we handle it in the rulebased solver. Our approach is less efficient computationally, but it ensures that all biochemical quantities are consistent and avoids having to make any assumptions about the relative speeds of processes.
Our approach allows us to model at a considerable level of detail. For example the conformation of NMDARs may be part of the biochemical model, allowing proteins in the PSD (e.g. calmodulin bound to calcium) to modulate the state of the channel [27]. We can also use one rulebased scheme to model both the presynapse and the postsynapse, which could help to understand transynaptic signalling via molecules such as endocannabinoids [4].
Our implementation of our algorithm is publicly available (KappaNEURON; http://github.com/davidcsterratt/KappaNEURON) and under development. The next major feature planned is making available to NEURON SpatialKappa’s capability of simulating rulebased models with voxelbased diffusion.
Appendix A Appendix
a.1 Kappa Simulation Method
To understand the asynchronous nature of the Kappa simulation method, we first illustrate Gillespie’s direct method [10] by applying it to the kinetic scheme description of a calcium pump shown in Fig. 3A. Here Ca represents intracellular calcium, P represents a pump molecule in the membrane, \ceP.Ca is the pump molecule bound by calcium and and are rate coefficients, which are rescaled to the variables and as explained in Section 2.5. To apply the Gillespie method to this scheme:

Compute the propensities of the reactions and

The total propensity is

Pick reaction with probability

Pick time to reaction , where is a random number drawn uniformly from the interval .

Goto 1
Kappa uses an analogous method, but applied to rules that are currently active. Both methods are eventbased rather than timestep based.
a.2 Justification for Throwing Away Events
To justify throwing away events occurring after a time step ending at , we need to show that the distribution of event times (measured from ) is the same in two cases:

The event time is drawn from an exponential distribution (for ), where is the propensity.

An event time is drawn as above. If , accept as the event time. If , throw away this event time and sample a new interval from an exponential distribution with a time constant of , i.e. . Set the event time to .
In the second case, the overall distribution is:
(18)  
Here we have used . Thus the distributions are the same in both cases.
Footnotes
References
 Bhalla, U.S., Iyengar, R.: Emergent properties of networks of biological signalling pathways. Science 283, 381–387 (1999)
 Cannon, R.C., O’Donnell, C., Nolan, M.F.: Stochastic ion channel gating in dendritic neurons: Morphology dependence and probabilistic synaptic activation of dendritic spikes. PLoS Comput. Biol. 68, e1000886 (2010)
 Carnevale, T., Hines, M.: The NEURON Book. Cambridge University Press, Cambridge, UK (2006)
 Castillo, P.E., Younts, T.J., Chávez, A.E., Hashimotodani, Y.: Endocannabinoid signaling and synaptic function. Neuron 761, 70–81 (2012)
 Chylek, L.A., Stites, E.C., Posner, R.G., Hlavacek, W.S.: Innovations of the rulebased modeling approach. In: Prokop, A., Csukás, B. (eds.) Systems Biology, pp. 273–300. SpringerVerlag (300)
 Collins, M.O., Husi, H., Yu, L., Brandon, J.M., Anderson, C.N.G., Blackstock, W.P., Choudhary, J.S., Grant, S.G.N.: Molecular characterization and comparison of the components and multiprotein complexes in the postsynaptic proteome. J. Neurochem. 97, 16–23 (2006)
 Danos, V., Feret, J., Fontana, W., Krivine, J.: Scalable simulation of cellular signaling networks. In: Shao, Z. (ed.) Programming Languages and Systems, Lecture Notes in Computer Science, vol. 4807, chap. 10, pp. 139–157. Springer, Berlin, Heidelberg (2007)
 Faas, G.C., Raghavachari, S., Lisman, J.E., Mody, I.: Calmodulin as a direct detector of Ca signals. Nat. Neurosci. 143, 301–304 (2011)
 Faeder, J., Blinov, M., Hlavacek, W.: Rulebased modeling of biochemical systems with BioNetGen. In: Maly, I.V. (ed.) Systems Biology, Methods in Molecular Biology, vol. 500, pp. 113–167. Humana Press (2009)
 Gillespie, D.: Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361 (1977)
 Harris, K.M., Kater, S.B.: Dendritic spines: Cellular specializations imparting both stability and flexibility to synaptic function. Annu. Rev. Neurosci. 17, 341–371 (1994)
 Hodgkin, A.L., Huxley, A.F.: A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. (Lond.) 117, 500–544 (1952)
 Kiehl, T.R., Mattheyses, R.M., Simmons, M.K.: Hybrid simulation of cellular behavior. Bioinformatics 203, 316–322 (2004)
 Lisman, J.E., Zhabotinsky, A.M.: A model of synaptic memory: a CaMKII/PP1 switch that potentiates transmission by organizing an AMPA receptor anchoring assembly. Neuron 312, 191–201 (2001)
 Markram, H., Lübke, J., Frotscher, M., Sakmann, B.: Regulation of synaptic efficiency by coincidence of postsynaptic APs and EPSPs. Science 275, 213–215 (1997)
 Martin, S.J., Grimwood, P.D., Morris, R.G.M.: Synaptic plasticity and memory: An evaluation of the hypothesis. Annu. Rev. Neurosci. 231, 649–711 (2000)
 Mattioni, M., Le Novère, N.: Integration of biochemical and electrical SignalingMultiscale model of the medium spiny neuron of the striatum. PLoS ONE 87, e66811 (2013)
 McDougal, R.A., Hines, M.L., Lytton, W.W.: Reactiondiffusion in the NEURON simulator. Front. Neuroinform. 728 (2013)
 Pepke, S., KinzerUrsem, T., Mihalas, S., Kennedy, M.B.: A dynamic model of interactions of Ca, calmodulin, and catalytic subunits of Ca/calmodulindependent protein kinase II. PLoS Comput. Biol. 62, e1000675 (2010)
 Pocklington, A.J., Cumiskey, M., Armstrong, J.D., Grant, S.G.N.: The proteomes of neurotransmitter receptor complexes form modular networks with distributed functionality underlying plasticity and behaviour. Mol. Syst. Biol. 21 (2006)
 Ray, S., Bhalla, U.S.: PyMOOSE: Interoperable scripting in python for MOOSE. Front. Neuroinform. 2 (2008)
 Smolen, P., Baxter, D.A., Byrne, J.H.: A model of the roles of essential kinases in the induction and expression of late longterm potentiation. Biophys. J. 908, 2760–2775 (2006)
 Sorokina, O., Sorokin, A., Armstrong, J.D.: Towards a quantitative model of the postsynaptic proteome. Mol. Biosyst. 7, 2813–2823 (2011)
 Sorokina, O., Sorokin, A., Armstrong, J.D., Danos, V.: A simulator for spatially extended kappa models. Bioinformatics pp. 3105–3106 (2013)
 Sterratt, D., Graham, B., Gillies, A., Willshaw, D.: Principles of Computational Modelling in Neuroscience. Cambridge University Press, Cambridge, UK (2011)
 Stiles, J.R., Bartol, T.M.: Monte Carlo methods for simulating realistic synaptic microphysiology using MCell. In: De Schutter, E. (ed.) Computational Neuroscience: Realistic Modeling for Experimentalists, chap. 4, pp. 87–127. CRC Press, Boca Raton, FL (2001)
 Urakubo, H., Honda, M., Froemke, R.C., Kuroda, S.: Requirement of an allosteric kinetics of NMDA receptors for spike timingdependent plasticity. J. Neurosci. 2813, 3310–3323 (2008)
 Zeng, S., Holmes, W.R.: The effect of noise on CaMKII activation in a dendritic spine during LTP induction. J. Neurophysiol. 1034, 1798–1808 (2010)