Stochastic lattice model of synaptic membrane protein domains

TL;DR

This paper develops a stochastic lattice model for synaptic membrane protein domains, revealing how stochastic reaction-diffusion dynamics explain molecular turnover, fluctuations, and heterogeneity observed in synaptic domains, beyond mean-field predictions.

Contribution

It introduces a detailed stochastic reaction-diffusion model that captures the collective behavior and fluctuations of synaptic membrane protein domains, extending previous mean-field approaches.

Findings

Quantitative agreement with mean-field models of nonlinear diffusion.

Discrepancies between stochastic and mean-field reaction dynamics.

Stochastic dynamics explain molecular turnover and spatial heterogeneity.

Abstract

Neurotransmitter receptor molecules, concentrated in synaptic membrane domains along with scaffolds and other kinds of proteins, are crucial for signal transmission across chemical synapses. In common with other membrane protein domains, synaptic domains are characterized by low protein copy numbers and protein crowding, with rapid stochastic turnover of individual molecules. We study here in detail a stochastic lattice model of the receptor-scaffold reaction-diffusion dynamics at synaptic domains that was found previously to capture, at the mean-field level, the self-assembly, stability, and characteristic size of synaptic domains observed in experiments. We show that our stochastic lattice model yields quantitative agreement with mean-field models of nonlinear diffusion in crowded membranes. Through a combination of analytic and numerical solutions of the master equation governing the reaction dynamics at synaptic domains, together with kinetic Monte Carlo simulations, we find substantial discrepancies between mean-field and stochastic models for the reaction dynamics at synaptic domains. Based on the reaction and diffusion properties of synaptic receptors and scaffolds suggested by previous experiments and mean-field calculations, we show that the stochastic reaction-diffusion dynamics of synaptic receptors and scaffolds provide a simple physical mechanism for collective fluctuations in synaptic domains, the molecular turnover observed at synaptic domains, key features of the observed single-molecule trajectories, and spatial heterogeneity in the effective rates at which receptors and scaffolds are recycled at the cell membrane. Our work sheds light on the physical mechanisms and principles linking the collective properties of membrane protein domains to the stochastic dynamics that rule their molecular~components.

Tables

Table 1. Table 1: Unless indicated otherwise, we use here the same reaction kinetics and values of the dimensionless rate constants as in Refs.  Haselwandter2011 ; Kahraman2016 , which correspond to model C in Ref.  Haselwandter2015 and are consistent with experiments on glycine receptors and gephyrin scaffolds Specht2008 ; Choquet2003 ; Triller2008 ; Triller2005 : ( m 1 , m 2 , β , μ ) = b ​ ( 0.4 , 10 , 0.5 , 0.7 ) subscript 𝑚 1 subscript 𝑚 2 𝛽 𝜇 𝑏 0.4 10 0.5 0.7 (m_{1},m_{2},\beta,\mu)=b(0.4,10,0.5,0.7) and ( r ¯ , s ¯ ) = ( 0.05 , 0.05 ) ¯ 𝑟 ¯ 𝑠 0.05 0.05 (\bar{r},\bar{s})=(0.05,0.05) , with the right column in the table showing the connection between the notation used here and in Refs.  Haselwandter2011 ; Haselwandter2015 . As in Ref.  Kahraman2016 , we fix the time units in our model by adjusting the rate of receptor endocytosis within the range of values estimated previously Haselwandter2011 ; Haselwandter2015 from experiments, which correspond to characteristic time scales ranging from seconds to hours Specht2008 ; Choquet2003 ; Triller2008 , to k 1 = 1 / 750 subscript 𝑘 1 1 750 k_{1}=1/750   s − 1 superscript s 1 \text{s}^{-1} so that our model reproduces the scaffold recovery time measured in FRAP experiments Choquet2003 ; Specht2008 ; Calamai2009 . The indicated rate constants enter our (stochastic and mean-field) reaction-diffusion model of synaptic domains through Eq. ( 6 ) with Eqs. ( 15 )–( 23 ).

Chemical reactions	Rate constants
$R\stackrel{k_1}{\to }{R}_b$	$k_1=b\approx$ $1.3\times {10}^{-3}$ s-1
$R_b\stackrel{k_2}{\to }R$	$k_2=m_1\frac{\overline{r}}{1-\overline{r}-\overline{s}}\approx 3.0\times {10}^{-5}$ s-1
$M_b+R\stackrel{k_3}{\to }{M}_b+R_b$	$k_3=\frac{m_1\overline{r}+m_2\overline{s}}{\overline{r}\left(1-\overline{r}-\overline{s}\right)}\approx 1.5\times {10}^{-2}$ s-1
$R_b+S\stackrel{k_4}{\to }R+S$	$k_4=b\frac{\overline{r}}{\overline{s}}\frac{1}{1-\overline{r}-\overline{s}}\approx 1.5\times {10}^{-3}$ s-1
$R_b+R+S\stackrel{k_5}{\to }2R+S$	$k_5=\frac{m_2}{\overline{r}}\frac{1}{1-\overline{r}-\overline{s}}\approx 3.0\times {10}^{-1}$ s-1
$S\stackrel{k_6}{\to }{S}_b$	$k_6=\beta \approx 6.7\times {10}^{-4}$ s-1
$S_b\stackrel{k_7}{\to }S$	$k_7=\beta \frac{\overline{s}}{1-\overline{r}-\overline{s}}\approx 3.7\times {10}^{-5}$ s-1
$M_b+S\stackrel{k_8}{\to }{M}_b+S_b$	$k_8=\frac{\mu }{1-\overline{r}-\overline{s}}\approx 1.0\times {10}^{-3}$ s-1
$S_b+2S\stackrel{k_9}{\to }3S$	$k_9=\frac{\mu }{\overline{s}}\frac{2}{1-\overline{r}-\overline{s}}\approx 4.1\times {10}^{-2}$ s-1

Equations

$$\frac{\partial P}{\partial t}=\sum_{{\bf m}}\!\,\bigl{[}W({\bf N}-{\bf m};{\bf m})P({\bf N}-{\bf m},t)-W({\bf N};{\bf m})P({\bf N},t)\bigr{]}\,,$$

$$W=W_{\text{react}}+W_{\text{diff}}$$

$$W_{\text{diff}}=W_{\text{diff}}^r+W_{\text{diff}}^s$$

$$W_{\text{diff}}^{(1,2;\alpha )}(\mathbf{N};\mathbf{m})$$

$$\delta \left(m_i+{ϵ}^{\alpha }\right)\delta \left(m_{i\pm 1}-{ϵ}^{\alpha }\right)\prod _{j\ne i,i\pm 1}\delta \left(m_j\right),$$

$$W_{\text{react}}=\sum _l{W}_{\text{react}}^{(l)},$$

$$W_{\text{react}}^{(l)}(\mathbf{N};\mathbf{m})=\sum _i{\mathcal{R}}_i^{(l)}\prod _{j\ne i}\delta \left(m_j\right),$$

$$\frac{d{R}_i}{dt}$$

$$\frac{d{S}_i}{dt}$$

$$⟨{\eta }_i^{(\alpha )}(t_1){\eta }_j^{(\alpha )}(t_2)⟩=K_{i,j}^{(\alpha ;2)}\delta (t_1-t_2),$$

$$\frac{\partial r}{\partial t}=F^r(r,s)-{\nu }_r\nabla \cdot {\mathbf{J}}^r,$$

$$\frac{\partial s}{\partial t}=F^s(r,s)-{\nu }_s\nabla \cdot {\mathbf{J}}^s,$$

$$\displaystyle Q_{i\pm n}(t)=\sum\limits_{k=0}^{\infty}\frac{\partial^{k}q}{\partial x^{k}}\bigg{|}_{x=ia}\frac{(\pm an)^{k}}{k!}\,,$$

$${\mathbf{J}}^r$$

$${\mathbf{J}}^s$$

$${\mathcal{R}}_i^{(1)}$$

$${\mathcal{R}}_i^{(2)}$$

$${\mathcal{R}}_i^{(6)}$$

$${\mathcal{R}}_i^{(7)}$$

$${\mathcal{R}}_i^{(3)}$$

$${\mathcal{R}}_i^{(8)}$$

$${\mathcal{R}}_i^{(4)}=$$

$${\mathcal{R}}_i^{(5)}=$$

$${\mathcal{R}}_i^{(9)}=$$

$$\frac{\partial r}{\partial t}$$

$$\frac{\partial s}{\partial t}$$

$$\frac{\partial (r+s)}{\partial t}=\nu {\nabla }^2(r+s),$$

$$\frac{dr}{dt}$$

$$\frac{ds}{dt}$$

$$\frac{dP(N^s,t)}{dt}=$$

$$+k_7\frac{1-N^s+ϵ}{ϵ}P(N^s-ϵ,t),$$

$$\frac{dP(N^s,t)}{dt}=$$

$$+{\bar{k}}_8\frac{1-N^s+ϵ}{ϵ}(N^s-ϵ)P(N^s-ϵ,t),$$

$$\frac{dP(N^s,t)}{dt}=$$

$$+k_9\frac{1-N^s+ϵ}{ϵ}\frac{(N^s-ϵ)(N^s-2ϵ)}{2}$$

$$\times P(N^s-ϵ,t).$$

$$N^s=0\stackrel{{\alpha }_0}{\to }{N}^s=ϵ\stackrel{{\alpha }_1}{\to }...\stackrel{{\alpha }_{C-1}}{\to }{N}^s=1,$$

Full text

I Introduction

The stability and plasticity of synapses are thought to play central roles in memory formation and learning McAllister2007 ; Citri2008 . In particular, neurotransmitter receptor molecules, concentrated in postsynaptic membrane domains along with scaffolds and other kinds of proteins Ziv2014 ; Salvatico2015 , are crucial for signal transmission across chemical synapses Citri2008 ; Legendre2001 ; Specht2008 ; Tyagarajan2014 . The strength of the transmitted signal depends on the number of receptors localized in synaptic domains Citri2008 ; Specht2008 , and regulation of the receptor number in synaptic domains provides one mechanism for postsynaptic plasticity Carroll2001 ; Shepherd2007 ; Kneussel2014 . With the advent of in vivo superresolution light microscopy Hell1994 ; Choquet2003 ; Huang2010 ; Kusumi2014 , the multiscale properties of synaptic domains—from the stochastic diffusion trajectories of individual synaptic receptors to the overall size and stability of synaptic domains—can now be studied quantitatively Choquet2013 ; Ziv2014 ; Kneussel2014 ; Salvatico2015 . A central discovery Choquet2013 ; Ziv2014 ; Kneussel2014 ; Salvatico2015 here is that synaptic receptors Choquet2003 ; Triller2008 ; Triller2005 , as well as their associated scaffolds Okabe1999 ; Gray2006 ; Calamai2009 ; Ziv2014 , turn over rapidly, with individual molecules entering and leaving synaptic domains on typical time scales as short as seconds. In contrast, synaptic domains of a well-defined characteristic size can persist over months or even longer periods of time Trachtenberg2002 ; Grutzendler2002 , which may Citri2008 ; Carroll2001 ; Shepherd2007 ; Specht2008 ; Choquet2013 ; Ziv2014 ; Kneussel2014 ; Salvatico2015 constitute part of the cellular basis for memory formation and learning.

To explain the stability and characteristic size of synaptic domains in the face of rapid molecular turnover, a number of phenomenological models for membrane domain formation Sekimoto2009 ; Holcman2006 ; Earnshaw2006 ; Burlakov2012 ; Czondor2012 ; Shouval2005 —based on, for instance, balancing receptor fluxes into and out of synaptic domains or spatially varying effective reaction and diffusion rates—have been proposed. Through an interplay between quantitative experiments on minimal model systems lacking most of the synaptic machinery (for instance, single transfected fibroblast cells) and theoretical modeling it has been demonstrated Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 that the reaction and diffusion properties of receptors and their associated scaffolds at the cell membrane are sufficient for the self-assembly of stable synaptic receptor-scaffold domains of the characteristic size observed in neurons. In particular, the presence of a presynaptic terminal is not essential for the self-assembly of stable synaptic receptor-scaffold domains. The reaction-diffusion model of synaptic domains Haselwandter2011 ; Haselwandter2015 describing these experiments explains, based on a reaction-diffusion (Turing) instability Turing1952 of the mean-field equations governing receptor-scaffold reaction-diffusion dynamics, how interactions between receptors and their associated scaffolds, together with the diffusion properties of each molecule species at the cell membrane, are sufficient for the spontaneous formation, stability, and characteristic size of synaptic domains. For molecular reaction and diffusion rates consistent with experimental measurements on synaptic receptors and scaffolds Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 , the reaction-diffusion model yields Haselwandter2011 ; Haselwandter2015 , starting from random initial conditions, the self-assembly, stability, and characteristic size of synaptic domains observed in neurons. Conversely, it has been shown Calamai2009 ; Haselwandter2011 ; Haselwandter2015 that self-assembly of synaptic domains can be prevented in both experiment and theory through, for instance, selective modification of the reaction properties of scaffolds.

In common with other membrane protein domains Lang2010 ; Simons2010 ; Rao2014 ; Recouvreux2016 , synaptic domains are characterized Choquet2003 ; Specht2008 ; Triller2005 ; Triller2008 ; Ribrault2011 ; Choquet2013 ; Ziv2014 ; Salvatico2015 by low protein copy numbers ($\approx 10$–1000) and protein crowding. Previous work Butler2009 ; Butler2011 ; Cao2014 ; Lugo2008 ; Erban2009 ; Hecht2010 ; Erban2007 ; Cohen2005 ; Wylie2006 ; Samoilov2006 ; McKane2004 ; McKane2005 ; Gillespie1976 ; Muratov2005 ; Muratov2007 ; Lizana2008 ; Fanelli2010 ; Fanelli2013 suggests that the coupling between molecular noise and the nonlinear reaction-diffusion dynamics induced by protein crowding can lead to a rich interplay between fluctuations and deterministic dynamics at synaptic domains. Indeed, experiments Ribrault2011 ; Choquet2013 and theoretical modeling Shouval2005 ; Holcman2006 ; Sekimoto2009 ; Burlakov2012 ; Czondor2012 indicate that synaptic domains undergo collective fluctuations that may affect synaptic signaling. In a previous paper Kahraman2016 we demonstrated that the stochastic lattice model associated with the mean-field reaction-diffusion dynamics at synaptic domains Haselwandter2011 ; Haselwandter2015 yields, for molecular reaction and diffusion rates consistent with experimental measurements on synaptic receptors and scaffolds Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 , emergence of synaptic domains in the presence of rapid stochastic turnover of individual molecules, provides a quantitative link between the molecular noise inherent in reaction-diffusion processes and collective fluctuations in synaptic domains, and allows prediction of the stochastic dynamics of individual synaptic receptors and scaffolds at the cell membrane. In the present article we build on this previous work Kahraman2016 to provide a detailed discussion of the stochastic lattice model of receptor-scaffold reaction-diffusion dynamics at synaptic domains Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 and its relation to the corresponding mean-field model. We show Kahraman2016 that molecular noise can yield substantial deviations from mean-field results for the receptor-scaffold reaction-diffusion dynamics at synaptic domains, and that stochastic lattice models can be employed successfully to provide quantitative insights into the single-molecule and collective dynamics of membrane protein domains Lang2010 ; Simons2010 ; Rao2014 ; Recouvreux2016 ; Choquet2003 .

This article is organized as follows. We first summarize, in Sec. II, the stochastic lattice model of receptor-scaffold reaction-diffusion processes at synaptic domains Kahraman2016 ; Haselwandter2015 , which is defined mathematically by a suitable master equation (ME), and its relation to the corresponding mean-field model Haselwandter2015 ; Haselwandter2011 formulated in accordance with the standard formalism of chemical dynamics Butler2011 ; Cao2014 ; Lugo2008 ; Erban2009 ; Butler2009 ; Hecht2010 ; Erban2007 ; Cohen2005 ; Wylie2006 ; Samoilov2006 ; McKane2004 ; McKane2005 ; Gillespie1976 ; Gillespie1977 ; Gillespie2013 ; Epstein1998 ; Walgraef1997 ; Cross2009 ; Cross1993 ; Murray2002 ; Meinhardt1982 ; Maini2001 . We then provide a detailed discussion of the relation between stochastic and mean-field results for the diffusion-only (see Sec. III) and reaction-only (see Sec. IV) systems. We derive analytic solutions of the ME for special cases of the reaction dynamics at synaptic domains, and carry out extensive kinetic Monte Carlo (KMC) simulations of the ME for the diffusion-only and reaction-only systems. Allowing for an interplay between reaction and diffusion processes at the cell membrane we explore, in Sec. V, collective fluctuations in synaptic domains Ribrault2011 ; Choquet2013 , the molecular turnover at synaptic domains measured in fluorescence recovery after photobleaching (FRAP) experiments Choquet2003 ; Specht2008 ; Calamai2009 , and the stochastic single-molecule dynamics at synaptic domains Choquet2003 ; Choquet2013 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Specht2008 ; Triller2005 ; Triller2008 . We conclude, in Sec. VI, with a summary and discussion of our key results. Appendices A and B provide mathematical details pertaining to our analytic solutions of the ME for the reaction-only system.

II Reaction-diffusion model of synaptic receptor-scaffold domains

In this section we summarize the stochastic lattice model of the reaction and diffusion dynamics of synaptic receptor ($R$) and scaffold ($S$) molecules developed in Refs. Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 , which we use throughout this article (see Fig. 1). In this model, the membrane is discretized into membrane patches (lattice sites). We assume that chemical reactions only take place among receptors or scaffolds occupying the same lattice site, with random hopping of receptors and scaffolds between nearest-neighbor lattice sites. We focus here on the most straightforward scenario of a 1D system of length $L$ with $K$ patches of size $a=L/K$. The 2D formulation of our model Haselwandter2011 ; Haselwandter2015 shows Kahraman2016 similar stochastic dynamics of synaptic domains as the 1D formulation we consider here. We denote the hopping rates of receptors and scaffolds at lattice site $i$ by $D_{i}^{r}/\tau_{r}$ and $D_{i}^{s}/\tau_{s}$, where $D_{i}^{r}(t)$ and $D_{i}^{s}(t)$ model spatiotemporal variations in the receptor and scaffold hopping rates.

Synaptic membrane domains are crowded with molecules Choquet2013 ; Ziv2014 , which is expected Triller2005 ; Specht2008 ; Triller2008 ; Ribrault2011 ; Salvatico2015 to affect diffusion and reaction processes at synaptic domains. To account for molecular crowding in our model, we impose Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 the constraint that the rates of all reaction and diffusion processes that increase the receptor or scaffold number at a lattice site $i$ are $\propto\left(1-N^{r}_{i}-N^{s}_{i}\right)$, where $N^{r}_{i}/\epsilon^{r}$ and $N^{s}_{i}/\epsilon^{s}$ are the occupation numbers of receptors and scaffolds at site $i$ with the normalization constants $\epsilon^{r}$ and $\epsilon^{s}$ so that, at each site, the number of receptors and scaffolds cannot increase beyond $1/\epsilon^{r}$ and $1/\epsilon^{s}$, respectively. As a result, we have $0\leqslant N^{r}_{i}+N^{s}_{i}\leqslant 1$ for all $i$. Analogous phenomenological models of crowding have been employed previously in a variety of different contexts Satulovsky1996 ; McKane2004 ; Lugo2008 ; Fanelli2010 ; Fanelli2013 . Based on recently-developed computational methodologies Flegg2012 ; Schoneberg2013 ; Gillespie2014 ; Johnson2014 ; Vijaykumar2015 ; Metzler2016 ; Jeon2016 for the description of reaction and diffusion processes at molecular scales, the simple model of crowding we focus on here could be connected to more detailed molecular models of the interactions between receptors and scaffolds.

In our stochastic reaction-diffusion model of synaptic domains Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 , the state of the system at time $t$ is completely characterized by the set of molecular occupation numbers $\mathbf{N}=\{\mathbf{N}^{\alpha}\}$ with $\alpha=r,s$, where $\mathbf{N}^{\alpha}(t)=(N^{\alpha}_{1}(t),N^{\alpha}_{2}(t),\cdots,N^{\alpha}_{K}(t))$. The stochastic dynamics of the system are governed by the ME Gardiner1985 ; Kampen1992

[TABLE]

where $P(\mathbf{N},t)$ is the probability that the system is in state $\mathbf{N}$ at time $t$ and $W(\mathbf{N};\mathbf{m})$ is the transition rate from state $\mathbf{N}$ to state $\mathbf{N}+\mathbf{m}$. Unless indicated otherwise, we use here Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 random initial conditions of $\mathbf{N}$ satisfying $0\leq N^{r}_{i}+N^{s}_{i}\leq 1$ for all $i$ with periodic boundary conditions.

The total transition rate $W$ in Eq. (1) can be written as

[TABLE]

in our reaction-diffusion model of synaptic domains Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 , where $W_{\text{react}}$ and $W_{\text{diff}}$ denote contributions to $W$ due to receptor and scaffold reaction and diffusion processes at the cell membrane. For the receptor and scaffold diffusion processes we have

[TABLE]

with $W_{\text{diff}}^{\alpha}=W_{\text{diff}}^{(1;\alpha)}+W_{\text{diff}}^{(2;\alpha)}$, in which the $W_{\text{diff}}^{(1,2;\alpha)}$ denote the receptor and scaffold transition rates for hopping from site $i$ to site $i\pm 1$,

[TABLE]

where the summation runs over the entire system, $\delta(x)$ is the Dirac-delta function, and the term $(1-N^{r}_{i\pm 1}-N^{s}_{i\pm 1})$ captures the effects of molecular crowding on receptor and scaffold diffusion Haselwandter2011 ; Haselwandter2015 . We use Dirac-delta functions, rather than Kronecker-delta functions, in Eq. (LABEL:eq:diffEx) in order to make the connection between the ME (1) and the corresponding mean-field equations more transparent, which amounts to replacing the summation in the ME (1) by an integral over all (continuous) ${\bf m}$ Haselwandter2007 . The factor of $1/\epsilon^{\alpha}$ in Eq. (LABEL:eq:diffEx) arises because we follow here the convention Kahraman2016 that $D_{i}^{\alpha}/\tau_{\alpha}$ is the hopping rate per molecule.

The contribution to $W$ due to reactions is given by

[TABLE]

in which each $W_{\text{react}}^{(l)}$ corresponds to a particular reaction among receptors or scaffolds. The $W_{\text{react}}^{(l)}$ take the general form

[TABLE]

where the summation runs over the entire system and, as in Eq. (LABEL:eq:diffEx), we use Dirac-delta functions so as to allow for continuous ${\bf m}$ in the ME (1). The $\mathcal{R}_{i}^{(l)}$ are dictated by the receptor or scaffold reaction dynamics Specht2008 ; Choquet2003 ; Triller2008 ; Triller2005 ; Haselwandter2011 ; Haselwandter2015 , and we return to their specific forms below.

To derive the mean-field equations associated with our stochastic lattice model Haselwandter2015 we introduce the continuum representations $R_{i}(t)$ and $S_{i}(t)$ of $N^{r}_{i}(t)$ and $N^{s}_{i}(t)$, respectively. Based on Eqs. (LABEL:eq:diffEx) and (6), the transition rates in our reaction-diffusion model can be directly extended to the continuous occupation numbers $R_{i}(t)$ and $S_{i}(t)$, allowing the ME (1) to be transformed Kampen1992 ; Haselwandter2007 ; Fox1991 ; Horsthemke1977 into the more tractable lattice Langevin equations Haselwandter2015

[TABLE]

where the $K_{i}^{(\alpha;1)}$ are the first moments of the contributions to $W$ changing the receptor or scaffold distribution in the system, and the Gaussian noise terms $\eta_{i}^{(\alpha)}$ have zero mean and covariance

[TABLE]

in which the $K_{i,j}^{(\alpha;2)}$ are the second moments of the contributions to $W$ changing the receptor or scaffold distribution in the system.

The continuum limit of the deterministic parts of the lattice Langevin equations (7) and (8) yields Haselwandter2011 ; Haselwandter2015 the mean-field equations

[TABLE]

with all parameters determined directly by the ME (1), where $r(x,t)$ and $s(x,t)$ are the continuum fields associated with $R_{i}(t)$ and $S_{i}(t)$ in Eqs. (7) and (8) with the noise terms set to zero,

[TABLE]

in which $Q_{i}\equiv R_{i},S_{i}$ and $q\equiv r(x,t),s(x,t)$, the polynomials $F^{\alpha}(r,s)$ in Eqs. (10) and (11) capture chemical reactions among receptors or scaffolds as in the standard formalism of chemical dynamics Butler2011 ; Cao2014 ; Lugo2008 ; Erban2009 ; Butler2009 ; Hecht2010 ; Erban2007 ; Cohen2005 ; Wylie2006 ; Samoilov2006 ; McKane2004 ; McKane2005 ; Gillespie1976 ; Gillespie1977 ; Gillespie2013 ; Epstein1998 ; Walgraef1997 ; Cross2009 ; Cross1993 ; Murray2002 ; Meinhardt1982 ; Maini2001 , the $\nu_{\alpha}=a^{2}/2\tau_{\alpha}$ are the receptor and scaffold diffusion coefficients, and the diffusion currents are given by

[TABLE]

where the $D^{\alpha}(x,t)$ denote the continuum representations Haselwandter2015 of $D^{\alpha}_{i}$ obtained via Eq. (12) with $Q_{i}\equiv D_{i}^{\alpha}$ and $q\equiv D^{\alpha}(x,t)$. The diffusion currents in Eqs. (LABEL:eq:dif_curR) and (LABEL:eq:dif_curS) follow directly Haselwandter2015 from the random hopping of receptors and scaffolds with rates $D^{\alpha}_{i}/\tau_{\alpha}$ together with the constraint that the rates of diffusion processes locally increasing the molecule occupancy at a given site $i$ are $\propto\left(1-N^{r}_{i}-N^{s}_{i}\right)$, as captured by Eq. (LABEL:eq:diffEx). Nonlinear crowding terms equivalent to those in Eqs. (LABEL:eq:dif_curR) and (LABEL:eq:dif_curS) have been studied previously in population biology Satulovsky1996 ; McKane2004 ; Lugo2008 and in the context of general models of non-Fickian diffusion Fanelli2010 ; Fanelli2013 . Membrane-bound reaction-diffusion systems similar to Eqs. (10) and (11) occur in a variety of different contexts Haselwandter2011 ; Haselwandter2015 ; Hecht2010 ; Levine2005 ; Marenduzzo2013 ; Vandin2016 ; Camley2017 .

We use glycine receptors and gephyrin scaffolds Ziv2014 ; Salvatico2015 ; Ribrault2011 ; Choquet2013 ; Specht2008 ; Triller2008 ; Triller2005 as a model system to fix the reaction kinetics and diffusion coefficients in our reaction-diffusion model of synaptic receptor-scaffold domains. We summarize here the pertinent reaction-diffusion dynamics, and refer the interested reader to Refs. Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 for a more detailed discussion of how the reaction-diffusion model considered here relates to the experimental phenomenology of glycine receptors and gephyrin. We first note that, at the lowest order, receptors and scaffolds may be randomly removed from the cell membrane via endocytosis as well as randomly inserted into the cell membrane, resulting in the reactions $R\xrightarrow{k_{1}}R_{b}$, $R_{b}\xrightarrow{k_{2}}R$, $S\xrightarrow{k_{6}}S_{b}$, and $S_{b}\xrightarrow{k_{7}}S$ (Fig. 1). In these expressions, $R$ and $S$ represent receptor and scaffold molecules at the membrane, while $R_{b}$ and $S_{b}$ stand for receptor and scaffold molecules in the cytoplasmic “bulk” of the cell Haselwandter2011 ; Haselwandter2015 , with the $k_{l}$ denoting rate constants. The resulting transition rates are given by Eq. (6) with Haselwandter2015 ; Kahraman2016

[TABLE]

yielding the additive contributions $-k_{1}r$ and $k_{2}(1-r-s)$ to $F^{r}$ in Eq. (10), and $-k_{6}s$ and $k_{7}(1-r-s)$ to $F^{s}$ in Eq. (11). The rate constants $k_{l}$ in Eqs. (15)–(18) [as well as the rate constants in Eqs. (19)–(23); see below] are scaled by $1/\epsilon^{\alpha}$ because we use the convention Kahraman2016 that $k_{l}$ denotes the rate of removal from/insertion into the cell membrane per molecule.

Furthermore, we note Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 that removal of receptors or scaffolds from the cell membrane may be facilitated by some mechanism that involves a temporary increase in the local crowding of the cell membrane, $M_{b}+R\xrightarrow{k_{3}}M_{b}+R_{b}$ and $M_{b}+S\xrightarrow{k_{8}}M_{b}+S_{b}$ (Fig. 1), where $M_{b}$ denotes an auxiliary bulk molecule. Because we account here for the effects of molecular crowding, these reactions yield contributions to the stochastic lattice model that are distinct from $R\xrightarrow{k_{1}}R_{b}$ and $S\xrightarrow{k_{6}}S_{b}$, and result in the terms Haselwandter2015 ; Kahraman2016

[TABLE]

in Eq. (6). Equations (19) and (20) imply the additive contributions $-k_{3}(1-r-s)r$ to $F^{r}$ in Eq. (10) and $-k_{8}(1-r-s)s$ to $F^{s}$ in Eq. (11), respectively.

Finally, we note that, as discussed previously Haselwandter2011 ; Haselwandter2015 , key experimental features of the reaction dynamics of glycine receptors and gephyrin for self-assembly of synaptic domains Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 are Specht2008 ; Choquet2003 ; Triller2008 ; Triller2005 ; Haselwandter2011 ; Haselwandter2015 that gephyrin can transiently bind glycine receptors as well as other gephyrin molecules, with experiments and theory suggesting Choquet2003 ; Calamai2009 ; Kahraman2016 ; Haselwandter2011 ; Haselwandter2015 that trimerization of gephyrin is a key reaction for self-assembly of synaptic domains. Allowing for the same order of receptor reactions as scaffold reactions, these considerations suggest Haselwandter2015 the reactions $R_{b}+S\xrightarrow{k_{4}}R+S$, $R_{b}+R+S\xrightarrow{k_{5}}2R+S$, and $S_{b}+2S\xrightarrow{k_{9}}3S$ (Fig. 1), resulting in the terms Haselwandter2015 ; Kahraman2016

[TABLE]

in Eq. (6). Equations (21) and (22) yield the additive contributions $k_{4}(1-r-s)s$ and $k_{5}(1-r-s)rs$ to $F^{r}$ in Eq. (10), and Eq. (23) implies the additive contribution $k_{9}(1-r-s)s^{2}/2$ to $F^{s}$ in Eq. (11).

Unless indicated otherwise, we use here the same values of $k_{l}$ as in Ref. Kahraman2016 (see Table 1) which, as discussed in Refs. Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 , are consistent with experiments on synaptic domains formed by glycine receptors and gephyrin scaffolds. Similarly, we use Kahraman2016 , unless indicated otherwise, the diffusion coefficients $\nu_{r}=10^{2}\nu_{s}=10^{-2}\mu\text{m}^{2}/\text{s}$, with the corresponding hopping rates $1/\tau_{\alpha}=2\nu_{\alpha}/a^{2}$ in Eq. (1), consistent with experiments on glycine receptors and gephyrin scaffolds Specht2008 ; Choquet2003 ; Triller2008 ; Triller2005 ; Calamai2009 ; Meier2001 ; Haselwandter2011 ; Haselwandter2015 . For simplicity, we set $\epsilon^{r}=\epsilon^{s}\equiv\epsilon$ and $D^{r}=D^{s}=1$ throughout this article. Distinct values of $\epsilon^{r}$ and $\epsilon^{s}$ could be used to provide a more detailed model of the receptor and scaffold numbers at synaptic domains, while spatiotemporal variations in $D^{r}$ or $D^{s}$ could be used Haselwandter2011 ; Haselwandter2015 , for instance, to model the effects of pre- and postsynaptic interactions on receptor or scaffold diffusion. Unless indicated otherwise, we set $\epsilon=1/100$ and $a\approx 80$ nm so that Erban2009 ; Cao2014 the membrane patch size is smaller than the expected typical size of synaptic domains Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 but large enough to accommodate multiple receptors and scaffolds, with size $\approx 5$–$10$ nm for glycine receptors and gephyrin Kim2006 ; Du2015 .

As discussed in Sec. IV, we obtained exact analytic solutions of the ME (1) for reaction-only systems involving a subset of the reactions in Eqs. (15)–(23). We supplemented these exact analytic solutions for general reaction schemes through direct numerical solutions of the ME (1), for which we used the Euler method. Furthermore, we carried out KMC simulations of the ME (1) Gillespie1976 ; Gillespie1977 ; Gillespie2013 for diffusion-only (see Sec. III), reaction-only (see Sec. IV), and reaction-diffusion (see Sec. V) systems employing the “spatial next reaction” method described in Ref. Elf2003 . In our implementation of the spatial next reaction method Elf2003 we used Gillespie’s “direct” method Gillespie1976 to choose, at each lattice site, which receptors or scaffolds undergo reaction or hopping processes. On this basis we were able to track individual receptors and scaffolds in our KMC simulations. Finally, we numerically solved the mean-field equations (10) and (11) using standard methods Mathematica with the initial conditions, boundary conditions, and parameter values employed for the ME (1).

III Protein diffusion in crowded membranes

In this section we focus on diffusion-only systems described by the ME (1) with $W_{\textrm{react}}=0$ and $D^{r}_{i}=D^{s}_{i}=1$. For such systems, Eqs. (10) and (11) imply the mean-field diffusion equations Satulovsky1996 ; McKane2004 ; Lugo2008 ; Fanelli2010 ; Fanelli2013 ; Haselwandter2011 ; Haselwandter2015

[TABLE]

The nonlinear terms in the mean-field equations (24) and (25) result from molecular crowding (steric exclusion), and impede diffusion into crowded membrane regions. In line with experiments and large-scale computer simulations of crowded membranes Metzler2016 ; Jeon2016 , the nonlinear diffusion terms in Eqs. (24) and (25) have been shown Fanelli2010 to result in mean-square displacement curves that bear signatures of anomalous diffusion. Below, we first consider the special case $\nu_{r}=\nu_{s}\equiv\nu$ in Eqs. (24) and (25), for which the total molecule concentration of receptors and scaffolds, $r+s$, obeys the standard (linear) diffusion equation

[TABLE]

and then discuss more complex scenarios corresponding to $\nu_{r}\neq\nu_{s}$.

III.1 Identical receptor and scaffold diffusion coefficients

In this section we focus on diffusion-only systems with the diffusion coefficients of receptors and scaffolds being equal to each other, $\nu_{r}=\nu_{s}\equiv\nu=0.01$ $\mu$m${}^{2}/$s (see Fig. 2). As described in Sec. II, we expect that $\nu_{r}>\nu_{s}$ for synaptic receptors and scaffolds Triller2008 ; Choquet2013 ; Choquet2003 ; Triller2005 ; Okabe1999 ; Gray2006 ; Calamai2009 ; Ziv2014 ; Kneussel2014 , but we consider here the case $\nu_{r}=\nu_{s}$ for completeness. As initial conditions we use adjacent step-profiles of receptors and scaffolds, with receptors and scaffolds changing from $N^{r}_{i}=0$ or $N^{s}_{i}=0$ ($r=0$ or $s=0$) to $N^{r}_{i}=1$ or $N^{s}_{i}=1$ ($r=1$ or $s=1$) [see the inset of Fig. 2(a)]. We find that, for small enough $\epsilon$, the mean-field equations (24) and (25) are in quantitative agreement with averages over KMC simulations of the underlying ME (1) with $W_{\textrm{react}}=0$ and $D^{r}_{i}=D^{s}_{i}=1$. For instance, Fig. 2(a) shows excellent agreement between mean-field equations and averages over KMC simulations for $\epsilon=1/40$. As the value $\epsilon$ is increased, the discreteness of the molecular diffusion processes becomes increasingly important, and the mean-field approach begins to yield inaccurate results for the average molecule concentrations [see Fig. 2(b)]. However, we find that even with $\epsilon=1/8$, which corresponds to a maximum molecule occupancy per lattice site of only eight molecules, the mean-field equations (24) and (25) capture the shape of the average diffusion profiles. As expected from Eqs. (24) and (25) with Eq. (26), the temporal evolution of the total molecule concentration $N^{r}_{i}+N^{s}_{i}$ ($r+s$) takes the form of a Gaussian profile corresponding to standard (Fickian) diffusion with diffusion coefficient $\nu$. In contrast, the individual receptor and scaffold concentration profiles do not take the form of Gaussian concentration profiles, with diffusion of receptor and scaffold populations being hindered by steric constraints.

To explore the effect of molecular crowding on the diffusion of individual molecules, we used KMC simulations to compute the mean-square displacement (MSD) curves of receptors located at different initial positions in Fig. 2 [see Fig. 3(a)]. We find that receptors initially located at the center of a crowded membrane region show effective hopping rates that are reduced substantially compared to the case of free diffusion, resulting in a reduced initial slope of the MSD. In contrast, receptors initially located near the boundary of a crowded membrane region can easily diffuse into membrane regions with few receptors or scaffolds, resulting in an initial slope of the MSD that is only reduced slightly compared to the case of free diffusion. As $t\to\infty$, all molecules have the same effective diffusion coefficient $\nu\left(1-\langle N_{i}^{r}+N_{i}^{s}\rangle\right)$ independent of their initial location, where $\langle N_{i}^{r}+N_{i}^{s}\rangle$ is the average number of molecules per lattice site in the system. Thus, molecular crowding can initially yield spatially heterogeneous MSD curves that bear a signature of the initial location of the molecule under consideration, and asymptotically results in a reduced effective diffusion coefficient, with the magnitude of the reduction in the effective diffusion coefficient governed by the value of $\langle N_{i}^{r}+N_{i}^{s}\rangle$ [see Fig. 3(b)].

III.2 Distinct receptor and scaffold diffusion coefficients

Synaptic receptors are thought to diffuse more rapidly than their associated scaffolds Triller2008 ; Choquet2013 ; Choquet2003 ; Triller2005 ; Okabe1999 ; Gray2006 ; Calamai2009 ; Ziv2014 ; Kneussel2014 , and we therefore focus here on the case $\nu_{r}>\nu_{s}$. Similarly as in Sec. III.1 we find Kahraman2016 that, provided $\epsilon\lessapprox 1/10$, the mean-field equations (24) and (25) describing a diffusion-only system composed of receptors and scaffolds with $\nu_{r}/\nu_{s}>1$ are in quantitative agreement with averages over KMC simulations of the underlying ME (1), independent of the particular value of $\nu_{r}/\nu_{s}$ considered (see Fig. 4). We thus find that, even in situations where steric constraints severely limit the maximum receptor or scaffold occupancy per membrane patch and the diffusion profiles deviate substantially from Fickian diffusion, the nonlinear mean-field diffusion equations (24) and (25) Satulovsky1996 ; McKane2004 ; Lugo2008 ; Fanelli2010 ; Fanelli2013 successfully predict the average concentration profiles implied by the underlying stochastic model, with all parameters in Eqs. (24) and (25) determined directly by the ME (1).

As $\nu_{r}/\nu_{s}$ is increased, molecular crowding produces increasingly complex receptor and scaffold concentration profiles (Fig. 4). In particular, using adjacent step-profiles of receptors and scaffolds as initial conditions, we find that the slowly-diffusing scaffolds act as a partially permeable (and dynamic) barrier to receptor diffusion that can, for large enough $\nu_{r}/\nu_{s}$, lead to non-monotonic receptor concentration profiles [see the right panel of Fig. 4(a)]. This can be understood by noting that, with the initial conditions used for Fig. 4, the slowly-varying scaffold concentration profiles show a pronounced maximum for a sustained period of time. Equation (24) indicates that membrane regions with $\nabla^{2}s<0$ locally depress the receptor concentration profile, because steric constraints make it unfavorable for receptors to diffuse into membrane regions that are crowded with scaffolds. As shown in the right panel of Fig. 4(a), a maximum in the scaffold concentration profile can therefore result in a local (transient) minimum in the receptor concentration profile and, hence, a multimodal receptor concentration profile. We also note that, with increasing $\nu_{r}/\nu_{s}$, the receptor profile becomes increasingly asymmetric, with an increasingly pronounced tail away from the scaffolds. As a result, the total molecule concentration profile $N^{r}_{i}+N^{s}_{i}$ ($r+s$) also becomes increasingly asymmetric [see Fig. 4(b)]. As expected from Eqs. (24) and (25), the total molecule concentration profile $N^{r}_{i}+N^{s}_{i}$ ($r+s$), as well as the receptor and scaffold concentration profiles, do not take the form of Gaussian diffusion profiles if $\nu_{r}\neq\nu_{s}$.

Starting from step-like initial concentration profiles of receptors and scaffolds, we used KMC simulations to compute the MSD curves of receptors [see Fig. 5(a)] and scaffolds [see Fig. 5(b)] located at left, center, and right positions in the initial molecule distributions of receptors and scaffolds. Consistent with Fig. 3(a), we find that receptors and scaffolds initially located at the boundaries of the crowded membrane region are least constrained by molecular crowding, and show the largest MSDs. Receptors initially located close to the center of the crowded membrane region show the smallest initial MSDs of all receptors, because their diffusion is hindered by receptors in one direction and by scaffolds in the other direction. In contrast, scaffolds close to the domain center show similar initial MSDs as scaffolds initially located at the center of the scaffold distribution. This can be understood by noting that, even though the diffusion of scaffolds initially located close to the center of the crowded membrane region is hindered by both receptors and scaffolds, the rapid diffusion of receptors soon allows scaffolds to diffuse into membrane regions that were initially fully occupied by receptors. Similarly as in Fig. 3, the asymptotic properties of the receptor and scaffold MSD curves in Fig. 5 are set by the receptor and scaffold diffusion constants scaled by $\left(1-\langle N_{i}^{r}+N_{i}^{s}\rangle\right)$, independent of the initial location of receptors or scaffolds. However, the crossover time from the initial to the asymptotic properties of the MSD curves is sensitive to the initial conditions used.

IV Protein reaction dynamics in crowded membranes

In this section we focus on reaction-only systems described by the ME (1) with $W_{\textrm{diff}}=0$, and the corresponding mean-field equations

[TABLE]

where, as described in Sec. II, the polynomials $F^{r}(r,s)$ and $F^{s}(r,s)$ are formulated in accordance with the standard formalism of chemical dynamics Butler2011 ; Cao2014 ; Lugo2008 ; Erban2009 ; Butler2009 ; Hecht2010 ; Erban2007 ; Cohen2005 ; Wylie2006 ; Samoilov2006 ; McKane2004 ; McKane2005 ; Gillespie1976 ; Gillespie1977 ; Gillespie2013 ; Epstein1998 ; Walgraef1997 ; Cross2009 ; Cross1993 ; Murray2002 ; Meinhardt1982 ; Maini2001 . For simplicity we omit, in this section, indices labeling lattice sites, and denote the receptor and scaffold numbers by $N^{r,s}$, and the corresponding mean-field receptor and scaffold concentrations by $r(t)$ and $s(t)$, respectively.

Previous studies Samoilov2006 ; Erban2009 ; Gillespie1976 ; Gillespie1977 ; Gillespie2013 ; McKane2004 ; Lugo2008 ; Kahraman2016 have shown that, for finite molecule numbers, the deterministic mean-field descriptions used in standard models of chemical dynamics Epstein1998 ; Walgraef1997 ; Cross2009 ; Cross1993 ; Murray2002 ; Meinhardt1982 ; Maini2001 can fail to capture the average dynamics, as well as steady states, of the underlying stochastic reaction processes. As illustrated below, upper limits on the protein copy number due to molecular crowding, absorbing (non-fluctuating) states Samoilov2006 , bistability Samoilov2006 , and amplification of noise through nonlinear chemical reactions Erban2009 provide specific physical mechanisms yielding disagreement between MEs and mean-field equations. In Secs. IV.1–IV.3 we consider reaction processes among receptors or scaffolds of increasing complexity, inspired by the reaction dynamics at synaptic domains (see Sec. II). We show that it is, at least in some special cases, practical to directly solve, either analytically or numerically, the MEs describing the chemical reaction dynamics at synaptic domains. We test the validity of these direct solutions of the MEs using KMC simulations, and provide systematic comparisons with solutions of the corresponding mean-field equations describing the reaction kinetics at synaptic domains.

IV.1 Single chemical reactions

In this section we focus on the single—linear and nonlinear—chemical reactions $S_{b}\xrightarrow{k_{7}}S$, $S_{b}+S\xrightarrow{\bar{k}_{8}}2S$, and $S_{b}+2S\xrightarrow{k_{9}}3S$ considered in Sec. II. We investigate each one of these reactions in turn. The reaction $S_{b}+S\xrightarrow{\bar{k}_{8}}2S$ is thereby implicit Haselwandter2015 in the reaction $M_{b}+S\xrightarrow{k_{8}}M_{b}+S_{b}$ discussed in Sec. II. Following Sec. II, the mean-field equations associated with these reactions are given by

[TABLE]

with the underlying MEs

[TABLE]

We take the initial scaffold concentration at $t=0$ to be given by $N^{s}_{0}=s_{0}$ for the MEs and the mean-field equations, which then uniquely specifies the solutions of Eqs. (29)–(31) and Eqs. (32)–(34), respectively. Note that Eqs. (32)–(34) exhibit an absorbing state at $N^{s}=1$ with $P=1$ for $N^{s}=1$ and $P=0$ otherwise, with Eqs. (33) and (34) also exhibiting an absorbing state at $N^{s}=0$, and a further absorbing state at $N^{s}=\epsilon$ in Eq. (34). All transition rates that would allow the system to exit an absorbing state are, by definition, equal to zero. Fluctuations are therefore completely suppressed in an absorbing state, independent of the value of $\epsilon$ considered. The steady states of the mean-field equations (29)–(31) are given by $s=0$ or $s=1$.

The MEs (32)–(34) can be solved analytically, as follows. We first note that Eqs. (32)–(34) correspond to chains of Markov processes that irreversibly transform a system with $N_{0}^{s}/\epsilon$ scaffold molecules into a system with $1/\epsilon$ scaffold molecules. The longest possible chain of reactions is obtained with $N_{0}^{s}=0$, for which

[TABLE]

where $C=1/\epsilon$ and the rates $\alpha_{i}$, in which $i={0,1,\dots,C-1}$, are obtained from the first terms on the right-hand sides of the MEs (32)–(34):

[TABLE]

The probability distribution of jump times $t$ between two consecutive states $i$ and $i+1$ of the Markov processes in Eq. (35) is given by the distribution of waiting times of a Poisson process with rate $\alpha_{i}$,

[TABLE]

Equation (39) with Eqs. (36)–(38) allows direct calculation of the probability distribution of jump times between arbitrary states in Eq. (35) (see Appendix A). However, to calculate the mean jump time between two arbitrary states in Eq. (35) it is sufficient to note from Eq. (39) that the mean jump time from state $i$ to state $i+1$ is given by $\langle t\rangle_{i\to i+1}=1/\alpha_{i}$, implying a mean time

[TABLE]

to reach a state $q$ with scaffold occupancy $N^{s}=q\epsilon$ starting from an initial state $p=N^{s}_{0}/\epsilon$ (see Appendix B).

We find that the mean jump times predicted by Eq. (40) with Eqs. (36)–(38) are in quantitative agreement with averages over KMC simulations of the corresponding MEs (32)–(34), for all values of $\epsilon$ considered here (see Fig. 6). To quantify the extent to which stochastic and mean-field results are in agreement with each other in Fig. 6 we calculate, for each time point available from the MEs (32)–(34) in Fig. 6, the difference between average stochastic and deterministic results, divide this difference by the corresponding mean-field result, and average the resultant percentage difference between stochastic and mean-field results over all time points in Fig. 6 associated with a given reaction and value of $\epsilon$. For $\epsilon\gtrapprox 1/5$ we find discrepancies $\gtrapprox 10\%$ between the mean-field equations (29)–(31) and the MEs (32)–(34). The disagreement between Eqs. (29)–(31) and Eqs. (32)–(34) becomes increasingly pronounced with increasing order of the reaction. For instance, for the linear reaction $S_{b}\xrightarrow{k_{7}}S$, $\epsilon\lessapprox 1/200$ gives a discrepancy $\lessapprox 0.5\%$ between the mean-field and master equations, while the reaction $S_{b}+2S\xrightarrow{k_{9}}3S$ requires $\epsilon\lessapprox 1/5000$ for mean-field and master equations to be in similarly good agreement. In particular, for the linear reaction $S_{b}\xrightarrow{k_{7}}S$ the ME (32) yields, for large enough $\epsilon$, a more rapid temporal evolution of the system than the corresponding mean-field equation (29) [see Fig. 6(a)]. This can be understood by noting that the ME allows for discrete jump processes with the system reaching the state $N^{s}=1$, on average, in a finite amount of time given by Eq. (40), while the mean-field solution only reaches $s=1$ as $t\to\infty$. The same argument holds true for the nonlinear reactions in Fig. 6 [see Figs. 6(b) and 6(c)]. However, for these nonlinear reactions we also find that, for a substantial portion of the trajectory of the system, the mean jump time in the MEs (33) and (34) is increased compared to the corresponding mean-field results implied by Eqs. (30) and (31). These results illustrate that, already for very simple reaction dynamics, the molecular noise inherent in scaffold (and receptor) reaction dynamics can have subtle effects on average system properties if steric effects prohibit maximum molecule occupancies greater than hundreds of molecules, as is often the case for proteins in cell membranes.

IV.2 Competing chemical reactions

In this section we study minimal reaction dynamics with competing chemical reactions increasing and decreasing the molecule number in the system. As model systems we use some of the key scaffold reactions at synaptic domains (see Sec. II). Reaction schemes with competing chemical reactions may, on the one hand, exhibit absorbing states in which all fluctuations are suppressed. On the other hand, competing chemical reactions may also yield fluctuating steady states of the ME with a mean that, for $\epsilon>0$, does not necessarily coincide with the steady state implied by the corresponding mean-field model. We first consider reaction schemes for which the mean-field model predicts steady state(s) with $0<s<1$, while the ME yields absorbing state(s) at $N^{s}=0$ or $N^{s}=1$. Consistent with previous work Samoilov2006 we find that absorbing states can yield breakdown of the mean-field approach. We then consider reaction schemes with competing chemical reactions for which the MEs exhibit fluctuating steady states and no absorbing states. We find quantitative agreement between mean-field predictions and averages over the underlying MEs for linear reaction schemes, independent of the value of $\epsilon$ considered. However, we also find that the mean-field approach can break down, for $\epsilon>0$, even for very simple nonlinear reaction schemes exhibiting fluctuating steady states.

IV.2.1 Absorbing states

Consider a system composed of scaffolds undergoing the reactions $S\xrightarrow{k_{6}}S_{b}$ and $S_{b}+S\xrightarrow{k_{9}}2S$ discussed in Sec. II. Following Sec. II, the corresponding mean-field equation is given by

[TABLE]

with the underlying ME

[TABLE]

Equation (41) yields, for $k_{9}>k_{6}$, a stable steady state at $s=1-k_{6}/k_{9}$ and an unstable steady state at $s=0$, while Eq. (42) implies an absorbing state at $N^{s}=0$ with $P=1$ for $N^{s}=0$ and $P=0$ otherwise.

Direct solution of Eq. (41) for $k_{9}>k_{6}$ confirms that, provided $s>0$ initially, the mean-field solution approaches the stable steady state $s=1-k_{6}/k_{9}$ as $t\to\infty$ [see Fig. 7(a)]. KMC simulations of Eq. (42) show that, initially, the average $N^{s}$ also tends to approach the stable steady state of the mean-field model [Fig. 7(a)]. Over time, however, fluctuations gradually carry the average of the stochastic system, with $\epsilon>0$, away from the steady state of the mean-field system and towards the absorbing state. This can be seen quite clearly by following individual stochastic trajectories of the system, which can get irreversibly locked into the absorbing state $N^{s}=0$ [see Fig. 7(b)]. Thus, the asymptotic behavior of the mean-field equation (41) is dominated by its stable steady state at $s=1-k_{6}/k_{9}$, while the asymptotic behavior of the underlying stochastic system is dominated by the absorbing state at $N^{s}=0$. As $\epsilon$ is increased, averages over the ME (42) approach the absorbing state increasingly rapidly, with the magnitude of the average slope increasing exponentially with decreasing maximum molecule occupancy (decreasing $1/\epsilon$) [see inset in Fig. 7(a)]. We find that, depending on the value of $\epsilon$ considered, it can take a substantial amount of time until the absorbing state is reached in our KMC simulations of the ME (42). For instance, with $\epsilon=1/400$ and $k_{6}$ and $k_{9}$ both being of the order of $1$ s*-1*, only $\approx 3\%$ of the 2000 KMC trajectories in Fig. 7(a) are absorbed by $N^{s}=0$ for $t\lessapprox 4000$ s.

As a second example of a system with competing chemical reactions exhibiting an absorbing state, consider the reactions $M_{b}+S\xrightarrow{k_{8}}M_{b}+S_{b}$ and $S_{b}+S\xrightarrow{\bar{k}_{8}}2S$ (see also Sec. IV.1). Following Sec. II, the corresponding mean-field equation is given by

[TABLE]

with the ME

[TABLE]

The ME (IV.2.1) has two absorbing states, at $N^{s}=0$ and $N^{s}=1$. For $k_{8}\neq\bar{k}_{8}$, the mean-field equation (43) exhibits steady states at $s=0$ and $s=1$, with the steady state $s=0$ ($s=1$) being unstable (stable) for $k_{8}<\bar{k}_{8}$, and vice versa for $k_{8}>\bar{k}_{8}$. For $k_{8}=\bar{k}_{8}$, the right-hand side of the mean-field equation (43) is identical to zero, with any initial condition $s(0)=s_{0}$ corresponding to a steady state of the system.

Comparing KMC simulations of the ME (IV.2.1) with solutions of the mean-field equation (43) we find that, for $k_{8}=\bar{k}_{8}$, the mean-field equation (43) is in quantitative agreement with averages over the ME (IV.2.1), independent of the initial conditions used [see the left panel of Fig. 8(a)]. However, the average scaffold concentration fails to capture the asymptotic properties of individual KMC trajectories which, for long enough times, are absorbed by either $N^{s}=0$ or $N^{s}=1$, and hence do not fluctuate about the average scaffold concentration [see the right panel of Fig. 8(a)]. For $k_{8}\neq\bar{k}_{8}$, we find that the mean-field equation (43) fails to capture averages over the ME (IV.2.1) with $\epsilon>0$ obtained via direct numerical solution of the ME (IV.2.1) and KMC simulations [see the left panel of Fig. 8(b)]. Similarly as for $k_{8}=\bar{k}_{8}$, the properties of individual stochastic trajectories of the system are dominated, at long times, by the absorbing states $N^{s}=0$ or $N^{s}=1$ rather than the average scaffold concentration [see the right panel of Fig. 8(b)]. We have confirmed the results of our KMC simulations and our numerical solutions of the ME (IV.2.1) through analytic solution of the ME (IV.2.1) at steady state. Our results illustrate Samoilov2006 how, for $\epsilon>0$, absorbing states can yield breakdown of the mean-field reaction dynamics at synaptic domains, and produce pronounced discrepancies between individual stochastic trajectories of the system and ensemble averages.

IV.2.2 Fluctuating steady states

In this section we focus on linear and nonlinear scaffold reaction dynamics with no absorbing states but (fluctuating) steady states. We first consider the linear reactions $S\xrightarrow{k_{6}}S_{b}$ and $S_{b}\xrightarrow{k_{7}}S$ discussed in Sec. II, yielding the mean-field equation

[TABLE]

with the stable steady state $s=k_{7}/\left(k_{6}+k_{7}\right)$, and the ME

[TABLE]

To analytically determine the steady-state probability distribution(s) of the ME (IV.2.2), $P_{\infty}(N^{s})$, we generalize the generating-function approach for the solution of MEs described in Ref. Gardiner1985 to allow for steric constraints. To this end, we define a generating function

[TABLE]

where $N^{s}=n\epsilon$ and we only include terms corresponding to $0\leq N^{s}\leq 1$ because, by definition, $P_{\infty}(N^{s})=0$ outside this range. The generating function in Eq. (47) satisfies the following identities:

[TABLE]

It then follows that

[TABLE]

where we have used that, by definition, $P_{\infty}(N^{s})=0$ for $N^{s}<0$. Setting the left-hand side of the ME (IV.2.2) equal to zero, multiplying the right-hand side by ${\hat{s}}^{n}$, summing all terms from $n=0$ to $n=1/\epsilon$, and employing the identities in Eqs. (48)–(50), we find that, in the steady state(s) of Eq. (IV.2.2), the generating function obeys

[TABLE]

Equation (51) has the unique solution

[TABLE]

satisfying the normalization constraint $G(1)=1$. Hence, Eq. (IV.2.2) admits only one steady-state probability distribution. To determine the form of this distribution we note from Eq. (47) that

[TABLE]

Equations (52) and (53) imply that the steady-state probability distribution associated with Eq. (IV.2.2) takes the form of a binomial distribution,

[TABLE]

where $n=N^{s}/\epsilon$. Equation (54) shows that the mean scaffold occupancy at steady state is given by $\langle N_{s}\rangle=k_{7}/\left(k_{6}+k_{7}\right)$, which is identical to the steady-state value of $s$ predicted by the mean-field equation (45) independent of the value of $\epsilon$ in Eq. (IV.2.2). As exemplified by Eq. (54), fluctuating steady states allow, in contrast to absorbing states, fluctuations about the mean molecule occupancy. As expected from the above considerations, solution of the mean-field equation (45), the average $N^{s}$ implied by Eq. (54), and averages over KMC simulations of the ME (IV.2.2) yield, even for $\epsilon=1$, identical results for the steady state of the system [see Fig. 9(a)]. We also find that the mean-field model yields quantitative agreement with KMC simulations of the ME (IV.2.2) for transient regimes of the system, independent of the value of $\epsilon$ considered [Fig. 9(a)].

The two most straightforward nonlinear versions of the reaction scheme considered above are $S+M_{b}\xrightarrow{k_{8}}S_{b}$ and $S_{b}\xrightarrow{k_{7}}S$, and $S\xrightarrow{k_{6}}S_{b}$ and $S_{b}+S\xrightarrow{\bar{k}_{8}}2S$ (see Sec. II). However, both of these reaction schemes exhibit absorbing states, which makes them unsuitable for the purposes of the present discussion. Instead, to explore the properties of nonlinear fluctuating steady states in a simple model system, we consider the reactions $2S\xrightarrow{k}2S_{b}$ and $S_{b}\xrightarrow{k_{7}}S$. The first of these reactions may not be relevant for synaptic domains Haselwandter2015 ; Haselwandter2011 , but may occur in other contexts Erban2007 ; Walgraef1997 ; Cross2009 ; Cross1993 ; Murray2002 . Following similar steps as in Sec. II, we find that the mean-field equation associated with this reaction scheme is given by

[TABLE]

yielding a stable steady state at

[TABLE]

with the second steady state of Eq. (55) lying outside the physically relevant range $0\leq s\leq 1$. The underlying ME is given by

[TABLE]

We ascertain the properties of the ME (IV.2.2) using direct numerical solutions of Eq. (IV.2.2) as well as KMC simulations. We find quantitative agreement between direct numerical solutions of the ME (IV.2.2) and KMC simulations for all values of $\epsilon$ considered here [see Fig. 9(b)]. For $\epsilon\gtrapprox 1/10$, however, the mean-field equation (55) fails to predict the average scaffold concentrations implied by the ME (IV.2.2) [Fig. 9(b)]. As $\epsilon$ is decreased, averages over the steady-state scaffold occupancies implied by the ME (IV.2.2) approach the steady-state value of $s$ in Eq. (56), with agreement between stochastic and mean-field results for $\epsilon\lessapprox 1/10$ [see Fig. 9(c)].

IV.3 Reaction kinetics at synaptic domains

In this section we consider the complete model of the reaction kinetics at synaptic domains described in Sec. II, which exhibits coupled, nonlinear reactions among receptors or scaffolds. We have shown previously Kahraman2016 that, for a maximum molecule occupancy per membrane patch $1/\epsilon\approx 100$, which is the regime of $\epsilon$ relevant for cell membranes (see also Sec. II), the mean-field equations (10) and (11) fail to capture the temporal evolution as well as steady-state values of the average receptor and scaffold concentrations implied by the ME (1), with the mean-field system approaching its steady state approximately one order of magnitude slower than the stochastic system. We have confirmed these conclusions using direct numerical solutions of the ME (1), which we find to be in quantitative agreement with KMC simulations of the ME (1) (see Fig. 10).

To further explore the stochastic reaction dynamics at synaptic domains, we follow individual stochastic trajectories of the system (see Fig. 11). Inspection of individual stochastic trajectories shows that, even in the steady state of the system, the receptor and scaffold occupancies can undergo large fluctuations, with the fluctuations in the receptor occupancy being particularly pronounced. To quantify the correlation between receptor and scaffold fluctuations, we calculate the receptor-scaffold correlation function in the steady-state of the stochastic system,

[TABLE]

in which we set $t_{\mathrm{min}}=1.0\times 10^{4}$ s and $t_{\mathrm{max}}=1.3\times 10^{5}$ s, $\langle N^{r,s}\rangle$ are the average receptor and scaffold occupancies in the time interval $[t_{\mathrm{min}},t_{\mathrm{max}}]$ implied by our KMC simulations, and we employ periodic boundary conditions when computing $N^{s}(t)$. We selected the value $t_{\mathrm{min}}=1.0\times 10^{4}$ s in Eq. (LABEL:eq:correlation) so that, for all the values of $\epsilon$ considered here, the average receptor and scaffold occupancies in Fig. 10 have reached their steady-state values at $t=t_{\mathrm{min}}$.

For all the values of $\epsilon$ considered here, we find that the time of maximum correlation in Eq. (LABEL:eq:correlation) occurs for $\bar{t}>0$ (see Fig. 12). This indicates that, consistent with the roles of receptors and scaffolds as “inhibitors” and “activators” of increased molecule occupancies Haselwandter2011 ; Haselwandter2015 , fluctuations increasing the scaffold occupancy trigger increased receptor occupancies. In turn, increased receptor occupancies tend to inhibit increased receptor as well as scaffold occupancies Haselwandter2011 ; Haselwandter2015 , with fluctuations decreasing the scaffold occupancy precipitating decreased receptor occupancies. Taken together, Figs. 11 and 12 thus suggest that small fluctuations in the scaffold occupancy can trigger large changes in the receptor occupancy, with the feedback between receptor and scaffold occupancies effectively amplifying receptor fluctuations. Figure 12 also suggests that the time of maximum correlation between fluctuations in the receptor and scaffold occupancies grows with decreasing $\epsilon$. Finally, we note that calculation of the power spectrum Butler2009 ; Butler2011 of the fluctuating receptor and scaffold occupancies obtained from our KMC simulations does not yield a characteristic (non-zero) frequency of receptor and scaffold fluctuations.

The KMC trajectories in Fig. 11 suggest that, for large enough $\epsilon$, the steady-state receptor occupancy is bistable, fluctuating between $N^{r}\gtrapprox 0.8$ and a value of $N^{r}$ close to zero. To further quantify the fluctuations in the receptor occupancy, we used our KMC data to estimate the steady-state probability distribution of the receptor occupancy, marginalized over the scaffold distribution, for different values of $\epsilon$ [see Fig. 13(a)]. Consistent with Fig. 11 we find that, for the value $\epsilon\approx 1/100$ relevant for synaptic domains Kahraman2016 , the receptor occupancy is bistable, with maxima at $N^{r}\approx 0$ and $N^{r}\approx 0.8$. As $\epsilon$ is decreased below $\epsilon\approx 1/300$, we find a unique maximum (mode) of the probability distribution. However, even for $\epsilon\approx 7.8\times 10^{-6}$, which would correspond to a “well-mixed” membrane compartment holding up to $\approx 1.28\times 10^{5}$ receptors or scaffolds, we find that both the averages [see Fig. 13(b)] and modes [see Fig. 13(c)] of the steady-state receptor and scaffold probability distributions do not coincide with the steady-state receptor and scaffold occupancies implied by the mean-field equations (10) and (11), with the disagreement being more pronounced for receptors than scaffolds.

The above results show that for, $\epsilon>0$, the ME (1) and the corresponding mean-field equations (10) and (11) can yield qualitatively and quantitatively different results for the reaction kinetics at synaptic domains discussed in Sec. II Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 . We further characterize the dependence of the solutions of the ME (1) on $\epsilon$ by calculating the occurrence probabilities of the different chemical reactions comprising the reaction dynamics at synaptic domains (see Table 1), as a function of $\epsilon$ (see Fig. 14). For the receptors we find that, as $\epsilon$ is decreased, the occurrence probabilities of the reactions $R\xrightarrow{k_{1}}R_{b}$ and $R_{b}+R+S\xrightarrow{k_{5}}2R+S$ tend to decrease compared to the occurrence probabilities of the other reactions changing the receptor number in the system. For the scaffolds we find that, as $\epsilon$ is decreased, the occurrence probability of the reaction $S\xrightarrow{k_{6}}S_{b}$ tends to decrease compared to the occurrence probabilities of the other reactions changing the scaffold number in the system. As further discussed in Sec. V, calculation of the occurrence probabilities of receptor and scaffold reactions across synaptic domains permits insights into possible physical mechanisms underlying spatially heterogeneous reaction dynamics at synaptic domains Kahraman2016 ; Czondor2012 ; Blanpied2002 ; Earnshaw2006 .

V Stochastic reaction-diffusion dynamics at synaptic domains

In this section we consider the full reaction-diffusion dynamics at synaptic domains (see Sec. II). We first study, and contrast, the collective properties of synaptic domains obtained from our stochastic and mean-field models of receptor-scaffold reaction-diffusion dynamics at synaptic domains. We then use our stochastic lattice model to compute the occurrence probabilities of receptor and scaffold reactions across synaptic domains. We find that the reaction-diffusion processes at synaptic domains give rise to characteristic locations of receptor and scaffold insertion and removal at synaptic domains. Next, we use our stochastic lattice model to study global receptor and scaffold turnover at synaptic domains Kahraman2016 , and make comparisons to the results of FRAP experiments Choquet2003 ; Specht2008 ; Calamai2009 . Finally, we show Kahraman2016 that our stochastic lattice model of the reaction-diffusion dynamics at synaptic domains yields single-molecule trajectories consistent with experimental observations Choquet2003 ; Choquet2013 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Specht2008 ; Triller2005 ; Triller2008 , and allows prediction of the time scale of receptor trafficking between membrane regions inside and outside synaptic domains.

V.1 Collective properties of synaptic domains

KMC simulations of the ME (1) show that the stochastic reaction-diffusion dynamics considered here (see Sec. II) yield Kahraman2016 , starting from random initial conditions, in-phase receptor and scaffold domains [see Fig. 15(a)]. Using the same values of the reaction and diffusion rates as in the ME (1), the mean-field equations (10) and (11) yield self-assembly of stable receptor-scaffold domains of a similar characteristic wavelength $\approx 8.5$ $\mu$m as found in our KMC simulations [Fig. 15(a)], which is also consistent with the linear stability analysis of the mean-field equations (10) and (11) Haselwandter2011 ; Haselwandter2015 . Similarly as in 2D systems, where the reaction and diffusion rates used here yield Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 synaptic domains of a similar characteristic size as observed in experiments Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 , the 1D patterns obtained from both the ME (1) and the mean-field equations (10) and (11) are irregular, with substantial variation in the spacing between individual domains. Furthermore, the ME (1) yields, for $\epsilon\approx 1/100$, substantial fluctuations in the size and location of synaptic domains, over a time scale of several hours. Consistent with the results on the reaction-only system in Sec. IV.3, we find that, for $\epsilon\approx 1/100$, domain formation proceeds more rapidly in the stochastic lattice model than in the mean-field model Kahraman2016 , by approximately one order of magnitude. Thus, the molecular noise inherent in receptor and scaffold reaction and diffusion processes accelerates the self-assembly of synaptic domains.

Consistent with our mean-field model Haselwandter2011 ; Haselwandter2015 , we find Kahraman2016 that the ME (1) yields scaffold profiles across synaptic domains that tend to be more narrow than receptor profiles [Fig. 15(a)]. Scaffold domains are therefore more sharply defined than receptor domains, and it is convenient Kahraman2016 to specify domain boundaries in our stochastic lattice model by placing a threshold on the scaffold occupancy. We first remove small-scale fluctuations in the scaffold occupancy using a Savitzky-Golay filter Savitzky1964 (order 5, frame size 25), and then apply a threshold $\bar{N}^{s}$ on the scaffold occupancy of membrane patches to automate detection of domain boundaries. We fix the value of $\bar{N}^{s}$ by examining how the average receptor and scaffold numbers per synaptic domain, as well as the average number of synaptic domains in the system, change as $\bar{N}^{s}$ is varied [see Fig. 15(b)]. For small enough $\bar{N}^{s}$, we obtain a single domain in the system encompassing all receptors and scaffolds, yielding the global maxima of the average receptor and scaffold numbers per domain, and the global minimum of the average domain number. As $\bar{N}^{s}$ is increased, the average domain number tends to increase rapidly, because of the many small, transient domains produced by the noise in the system. Concurrently, we find a drop in the average receptor and scaffold numbers per synaptic domain. The average domain number peaks for a value of $\bar{N}^{s}$ yielding a maximum number of small, transient domains, at which point we also find local minima in the average receptor and scaffold numbers per synaptic domain. As $\bar{N}^{s}$ is increased further, small, transient domains are increasingly filtered out, producing a drop in the average domain number and an increase in the average receptor and scaffold numbers per synaptic domain. We find that the average scaffold number per synaptic domain peaks at $\bar{N}^{s}\approx 0.08$. Beyond this local maximum, the average scaffold number per synaptic domain only decreases gradually with increasing $\bar{N}^{s}$. At around $\bar{N}^{s}\approx 0.08$ we also find a local maximum in the average receptor number per synaptic domain, as well as a marked decrease in the rate at which the average domain number decreases with increasing $\bar{N}^{s}$ [Fig. 15(b)]. These results suggest that $\bar{N}^{s}=0.08$ provides a suitable threshold for quantifying the boundaries of the stochastic reaction-diffusion patterns implied by our KMC simulations of the ME (1), which is also confirmed by direct inspection of our KMC data [Fig. 15(a)].

Using $\bar{N}^{s}=0.08$ we find that the mean values and standard deviations of the receptor and scaffold numbers per synaptic domain are given by $745\pm 355$ and $226\pm 108$ for the KMC data in Fig. 15(a), respectively, with pronounced fluctuations in the in-domain receptor and scaffold population numbers, over a time scale of hours [see Fig. 16(a)]. In analogy to the receptor-scaffold correlation function defined for the reaction-only system in Eq. (LABEL:eq:correlation), we quantify the fluctuations in Fig. 16(a) by calculating the correlation function of the receptor and scaffold populations in a given synaptic domain,

[TABLE]

where the $N^{D}_{r,s}(t)$ denote the in-domain receptor and scaffold numbers at time $t$. The averages $\langle N^{D}_{r,s}\rangle$ are evaluated over the time window $t_{\mathrm{min}}\leq t\leq t_{\mathrm{max}}$. We choose $t_{\mathrm{min}}$ and $t_{\mathrm{max}}$ to correspond to the lifetime of a given synaptic domain at the membrane. In marked contrast to the results for the reaction-only system in Fig. 12, which yield a time of maximum correlation between fluctuations in the receptor and scaffold occupancies of the order of minutes, we find that the global maximum of $C^{D}_{r,s}(\bar{t})$ occurs at $\left|\bar{t}\right|<1$ s [see Fig. 16(b)]. Thus, our KMC simulations of the ME (1) suggest that diffusion strongly diminishes the time of maximum correlation between fluctuations in the receptor and scaffold populations at synaptic domains.

We now turn our focus to the occurrence probabilities of receptor and scaffold reactions across synaptic domains. In particular, we consider the synaptic domain delineated by black boundaries in Fig. 15(a), and first average the corresponding receptor and scaffold profiles over the lifetime of this domain [see Fig. 17(a)]. We compute, across the synaptic domain, the densities of all reactions inserting (removing) receptors and scaffolds into (from) the membrane, and average these densities over the lifetime of the domain. We denote the resultant reaction densities by $\phi^{+}_{R,S}$ ($\phi^{-}_{R,S}$), which we normalize by the total number of receptor and scaffold reactions that occurred for the spatial region and time interval considered here [see Fig. 17(b)]. We find that $\phi^{+}_{R}$ and $\phi^{\pm}_{S}$ trace the approximate domain profile with the maxima in $\phi^{+}_{R}$ and $\phi^{\pm}_{S}$ occurring close to the domain center, which also shows the largest concentration of receptors and scaffolds [Fig. 17(a)]. In contrast, $\phi^{-}_{R}$ shows a local minimum at the domain center, with local maxima of $\phi^{-}_{R}$ occurring just outside the synaptic domain. This model prediction is consistent with experimental observations Blanpied2002 suggesting that receptors are removed from the membrane in membrane regions adjacent to synaptic domains.

To further investigate the origin of the aforementioned qualitative differences in the profiles of $\phi^{+}_{R}$ and $\phi^{\pm}_{S}$, and $\phi^{-}_{R}$, across synaptic domains, we consider the occurrence probabilities of each individual receptor and scaffold reaction across synaptic domains [see Fig. 17(c)]. We find that the two dominant reactions for the receptors are $R_{b}+R+S\rightarrow 2R+S$ and $M_{b}+R\rightarrow M_{b}+R_{b}$, while $S_{b}+2S\rightarrow 3S$ and $S\rightarrow S_{b}$ are dominant for the scaffolds. In particular, $\phi^{-}_{R}$ is mainly set by the reaction $M_{b}+R\rightarrow M_{b}+R_{b}$, which shows a similar profile as $\phi^{-}_{R}$. Steric constraints reduce the effective rate of the reaction $M_{b}+R\rightarrow M_{b}+R_{b}$ in crowded membrane regions, yielding a local minimum of the occurrence probability of $M_{b}+R\rightarrow M_{b}+R_{b}$ at the domain center. However, receptors diffuse rapidly, and can therefore readily leave synaptic domains. As a result, a substantial number of receptors are removed via the reaction $M_{b}+R\rightarrow M_{b}+R_{b}$ in the (less crowded) membrane regions adjacent to synaptic domains, producing local maxima of the occurrence probability of $M_{b}+R\rightarrow M_{b}+R_{b}$ and, hence, $\phi^{-}_{R}$ close to the domain boundaries. Note that the reaction $M_{b}+S\rightarrow M_{b}+S_{b}$, which is the scaffold reaction analogous to $M_{b}+R\rightarrow M_{b}+R_{b}$, does not produce pronounced local maxima of $\phi^{-}_{S}$ close to the domain boundaries. This can be understood by noting that scaffolds diffuse less rapidly than receptors, and are therefore less likely to diffuse into less crowded membrane regions with a large effective rate of $M_{b}+S\rightarrow M_{b}+S_{b}$.

V.2 Molecular turnover

Synaptic domains have been observed Choquet2013 ; Ziv2014 ; Kneussel2014 ; Salvatico2015 ; Choquet2003 ; Triller2008 ; Triller2005 ; Okabe1999 ; Gray2006 ; Calamai2009 to be in a dynamic steady state, with rapid turnover of individual receptors and scaffolds. To study, in our stochastic lattice model, receptor and scaffold turnover at synaptic domains we proceed Kahraman2016 as in FRAP experiments Choquet2003 ; Specht2008 ; Calamai2009 , and label all receptors and scaffolds inside a synaptic domain at a given time. The fraction of unlabeled receptors and scaffolds inside the synaptic domain as a function of time, monitored starting from the time when the receptors and scaffolds inside the synaptic domain were labeled, then provides a measure of global molecular turnover at synaptic domains (see Fig. 18). The time scale in our reaction-diffusion model is set Haselwandter2011 ; Haselwandter2015 by the rate of receptor endocytosis. We adjust Kahraman2016 the rate of receptor endocytosis within the range of values suggested by experiments Haselwandter2011 ; Haselwandter2015 ; Choquet2003 ; Specht2008 ; Triller2008 to $k_{1}=1/750$ $\text{s}^{-1}$ (see Table 1) so that, consistent with FRAP experiments Choquet2003 ; Specht2008 ; Calamai2009 , $\approx 30\%$ of scaffolds, but $>95\%$ of receptors, are replaced, on average, within $5$ min. We find that, on average, $>99\%$ of receptors (scaffolds) are replaced within $\approx 7$ min ($\approx 80$ min). Single realizations of receptor and scaffold turnover at synaptic domains Kahraman2016 trace closely, in our model, the corresponding results obtained by averaging over several realizations (Fig. 18).

V.3 Single-molecule dynamics

In agreement with experiments Choquet2003 ; Specht2008 ; Triller2005 ; Triller2008 ; Ribrault2011 ; Choquet2013 ; Ziv2014 ; Salvatico2015 , our reaction-diffusion model predicts Kahraman2016 that receptors initially localized in synaptic domains tend to leave synaptic domains via diffusion, while scaffolds typically stay localized within synaptic domains over their lifetime at the membrane. Following the trajectories of individual receptors from insertion into the membrane until removal from the membrane, we find that receptors tend to exchange, via diffusion, between the inside and outside of synaptic domains over their lifetime at the membrane [see Fig. 19(a)]. In particular, our KMC simulations of the ME (1) imply that receptors traffick along the membrane between the inside and outside of synaptic domains over a time scale of seconds to minutes [see Fig. 19(b)], which is consistent with the switching times between intra- and extra-synaptic membrane regions typically found in single-molecule experiments on receptor diffusion at synaptic domains Choquet2003 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Triller2005 ; Specht2008 ; Triller2008 ; Choquet2013 .

VI Summary and conclusions

In neurons, neurotransmitter receptors are concentrated in membrane regions associated with postsynaptic domains McAllister2007 ; Ziv2014 ; Salvatico2015 ; Citri2008 ; Legendre2001 ; Specht2008 ; Tyagarajan2014 , which are enormously complex molecular assemblies. Experiments on minimal systems devoid of the molecular machinery commonly associated with postsynaptic domains McAllister2007 ; Citri2008 —such as single transfected fibroblast cells—have shown Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 that receptor-scaffold interactions, together with the diffusion properties of each molecule species at the cell membrane, are sufficient for the formation, stability, and characteristic size of synaptic domains observed in neurons. The observed self-assembly of synaptic domains of a stable characteristic size can be understood Haselwandter2011 ; Haselwandter2015 ; Turing1952 based on mean-field models of receptor-scaffold reaction-diffusion processes at the cell membrane. However, mean-field models cannot capture how the rapid stochastic dynamics of individual synaptic receptors and scaffolds Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Triller2005 ; Specht2008 ; Triller2008 relate Ribrault2011 ; Choquet2013 ; Ziv2014 ; Salvatico2015 to the observed collective properties of synaptic domains. In particular, experiments Ribrault2011 ; Choquet2013 and theoretical modeling Shouval2005 ; Holcman2006 ; Sekimoto2009 ; Burlakov2012 ; Czondor2012 suggest that synaptic domains undergo collective fluctuations that may affect synaptic signaling.

In this article we have provided a detailed discussion of a stochastic lattice model of receptor-scaffold reaction-diffusion processes at the cell membrane Haselwandter2011 ; Haselwandter2015 that yields Kahraman2016 emergence of synaptic domains in the presence of rapid stochastic turnover of individual molecules. Based on the reaction and diffusion properties of synaptic receptors and scaffolds suggested by previous experiments and mean-field calculations Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 , we have shown previously Kahraman2016 that this stochastic lattice model provides a simple physical mechanism for collective fluctuations in synaptic domains Ribrault2011 ; Choquet2013 , the molecular turnover observed at synaptic domains Choquet2003 ; Specht2008 ; Calamai2009 , key features of the observed single-molecule trajectories Choquet2003 ; Choquet2013 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Specht2008 ; Triller2005 ; Triller2008 , and spatially inhomogeneous receptor and scaffold lifetimes at the cell membrane Czondor2012 ; Blanpied2002 ; Earnshaw2006 . We have confirmed here these conclusions Kahraman2016 , and expanded and elaborated upon our previous results, using a combination of KMC simulations of our stochastic lattice model of receptor-scaffold reaction-diffusion processes at synaptic domains, and analytic and numerical solutions of the ME governing our stochastic lattice model.

For the diffusion-only system we find that, even for $1/\epsilon=8$, which corresponds to a maximum molecule occupancy per lattice site of only eight molecules, the average results of KMC simulations are in quantitative agreement with the corresponding mean-field model Haselwandter2011 ; Haselwandter2015 ; Satulovsky1996 ; McKane2004 ; Lugo2008 ; Fanelli2010 ; Fanelli2013 of nonlinear diffusion in crowded membranes. A possible origin for this agreement between the stochastic lattice model and the mean-field model, even for large $\epsilon$, is that, in the diffusion-only system, the molecule number is conserved, which constrains the fluctuations in the stochastic system Haselwandter2002 . The ME (1) with $W_{\textrm{react}}=0$ and the corresponding mean-field equations (10) and (11) with $F^{r}=F^{s}=0$ Haselwandter2011 ; Haselwandter2015 ; Satulovsky1996 ; McKane2004 ; Lugo2008 ; Fanelli2010 ; Fanelli2013 produce non-Gaussian and, in some cases, even non-monotonic diffusion profiles. In particular, compared to Fickian diffusion, crowding tends to yield less disperse molecule distributions Kahraman2016 , with the slowly-diffusing scaffolds acting as an effective barrier to the dispersal of the more rapidly diffusing receptors Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Triller2005 ; Specht2008 ; Triller2008 ; Choquet2013 ; Ziv2014 ; Salvatico2015 . Calculation of the MSD of individual receptors and scaffolds shows that, at finite times, the MSD bears a signature of the initial distribution of the receptors and scaffolds in the system. As the system approaches its steady state, with a homogeneous distribution of receptors and scaffolds, we recover the characteristics of the MSD associated with standard Fickian diffusion, provided that the diffusion coefficient is rescaled by a factor accounting for the average molecular crowding in the system.

As discussed in Sec. II, for synaptic domains formed by glycine receptors and gephyrin, a physically reasonable choice for the value of $\epsilon$ is $\epsilon\approx 1/100$ Kahraman2016 so that, for a size $a_{P}\approx 5$–$10$ nm of glycine receptors and gephyrin Kim2006 ; Du2015 , the membrane patch size $a\approx 80$ nm is smaller than the expected typical size of synaptic domains Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 . As $\epsilon\to 0$, one generally expects that averages over the ME (1) coincide with the solutions of the corresponding mean-field equations (10) and (11) Haselwandter2011 ; Haselwandter2015 . As summarized above, agreement between stochastic lattice model and mean-field model is already obtained, for the diffusion-only system, with $\epsilon\approx 1/10$. In contrast, for the reaction-only system we find that, for physically relevant values of $\epsilon$, averages over the stochastic lattice model do, in general, not coincide with solutions of the corresponding mean-field equations. Exceptions to this conclusion are provided by purely linear reaction schemes with only a single linear reaction or competing linear reactions yielding a fluctuating steady state of the system (and no absorbing state). Indeed, our exact analytic solutions of the ME show that, in the case of competing linear reactions with a fluctuating steady state, agreement between ME and mean-field equations can even be obtained with $\epsilon=1$, which corresponds to a maximum molecule occupancy per membrane patch of only one molecule.

Consistent with previous theoretical studies of stochastic reaction-diffusion systems Samoilov2006 ; Erban2009 ; Gillespie1976 ; Gillespie1977 ; Gillespie2013 ; McKane2004 ; Lugo2008 ; Kahraman2016 , we find that molecular noise generally plays a central role in the reaction dynamics at synaptic domains. In particular, for the nonlinear reactions thought to be relevant for synaptic domains (see Sec. II) we find, for physically relevant values of $\epsilon$, disagreement between solutions of the ME (1) with $W_{\textrm{diff}}=0$ and the corresponding mean-field equations (10) and (11) with $\nu_{r}=\nu_{s}=0$ Haselwandter2011 ; Haselwandter2015 ; Satulovsky1996 ; McKane2004 ; Lugo2008 ; Fanelli2010 ; Fanelli2013 . For instance, trimerization of scaffolds, $S_{b}+2S\rightarrow 3S$, is found Choquet2003 ; Calamai2009 ; Kahraman2016 ; Haselwandter2011 ; Haselwandter2015 to be crucial for the self-assembly of stable synaptic domains composed of glycine receptors and gephyrin. Even in a very simple reaction scheme, in which only this single reaction is considered, a value $\epsilon\lessapprox 1/5000$ is required to produce reasonable agreement between the stochastic model and the mean-field model. Such a value of $\epsilon$ implies a membrane patch size $a\approx 71a_{P}\approx 350$–$710$ nm for glycine receptors and gephyrin Kim2006 ; Du2015 . The assumption of a well-mixed system is not expected to be warranted over such large membrane patch sizes, resulting in breakdown of the mean-field approach Epstein1998 ; Walgraef1997 ; Cross2009 ; Cross1993 ; Murray2002 ; Meinhardt1982 ; Maini2001 even for single (nonlinear) chemical reactions. In addition to the amplification of noise through nonlinear chemical reactions with steric constraints Kahraman2016 ; Erban2009 , illustrated by $S_{b}+2S\rightarrow 3S$, we also Samoilov2006 find that absorbing (non-fluctuating) states and bistability provide physical mechanisms yielding disagreement between stochastic and mean-field models of the reaction dynamics at synaptic domains. In the case of competing absorbing states we find that the mean-field model can produce quantitative agreement with averages over the stochastic model while failing to reproduce the asymptotic properties of individual realizations of the stochastic system, with distinct stochastic trajectories of the system being trapped in distinct absorbing states.

Comparing stochastic and mean-field results for the complete reaction dynamics at synaptic domains (see Sec. II) we find Kahraman2016 that, for $\epsilon\approx 1/100$, the mean-field equations (10) and (11) with $\nu_{r}=\nu_{s}=0$ fail to capture the temporal evolution as well as steady-state values of the average receptor and scaffold occupancies implied by the ME (1) with $W_{\textrm{diff}}=0$, with the average stochastic dynamics being approximately one order of magnitude faster than the mean-field dynamics. Calculation of the marginal steady-state probability distribution of the receptor occupancy shows that one origin for this discrepancy between stochastic and mean-field solutions lies, for $\epsilon\gtrapprox 1/300$, in the bistability of the stochastic system. For $\epsilon\lessapprox 1/300$, we find a single mode of the marginal steady-state probability distribution of the receptor occupancy. However, our KMC simulations show that, even for $\epsilon\approx 7.8\times 10^{-6}$, the averages as well as the modes of the receptor and scaffold occupancies can disagree substantially with mean-field results. We find that, in our stochastic model, the receptor and scaffold occupancies can both undergo large fluctuations, with the fluctuations in the receptor occupancy being particularly pronounced. Calculation of the receptor-scaffold correlation function shows that, consistent with the roles of receptors and scaffolds as inhibitors and activators of increased molecule concentrations Haselwandter2011 ; Haselwandter2015 , fluctuations increasing (decreasing) the scaffold occupancy tend to lead to an increase (a decrease) in the receptor occupancy, with the time of maximum correlation between receptor and scaffold fluctuations being of the order of minutes. We find that the time of maximum correlation between receptor and scaffold fluctuations tends to increase with decreasing $\epsilon$.

Finally, we considered in this article the receptor-scaffold reaction-diffusion dynamics at synaptic domains suggested by previous experiments and mean-field calculations Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 . We find that, starting from random initial conditions, the stochastic lattice model described by the ME (1) yields Kahraman2016 spontaneous formation of in-phase receptor and scaffold domains with a characteristic wavelength consistent with the corresponding mean-field results implied by Eqs. (10) and (11). We have shown previously Haselwandter2011 ; Haselwandter2015 ; Kahraman2016 that, in 2D, the mean-field equations (10) and (11) and the corresponding ME (1) yield, for the reaction and diffusion rates suggested by experiments Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 and used here (see Sec. II), synaptic domains of a similar characteristic size as observed in experiments on neurons and transfected fibroblast cells Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 . Consistent with our results for the reaction-only system, we find Kahraman2016 , for $\epsilon\approx 1/100$, that molecular noise accelerates synaptic domain formation by approximately one order of magnitude compared to mean-field dynamics while producing, over a time scale of several hours, substantial fluctuations in the size and location of synaptic domains. These results illustrate the potential importance of stochastic effects Choquet2013 when describing synaptic domains. We also find that, consistent with experimental observations Ribrault2011 ; Choquet2013 , the molecular noise induced by the underlying reaction and diffusion dynamics of synaptic receptors and scaffolds can produce Kahraman2016 collective fluctuations in synaptic domains. Calculation of the correlation function of the in-domain receptor and scaffold population numbers shows that, in contrast to our results for the reaction-only system, the reaction-diffusion system only yields a very short time of maximum correlation between receptor and scaffold fluctuations. Thus, our results suggest that diffusion strongly diminishes the correlation time between fluctuations in the receptor and scaffold populations at synaptic domains Ribrault2011 ; Choquet2013 .

We find Kahraman2016 that, in both our mean-field and stochastic models, scaffold profiles across synaptic domains tend to be more narrow than receptor profiles. Scaffold domains therefore tend to be more sharply defined than receptor domains, and we quantify Kahraman2016 domain boundaries in our stochastic lattice model by placing a threshold on the scaffold occupancy per membrane patch. Our stochastic lattice model allows us to compute the occurrence probabilities of receptor and scaffold reactions across synaptic domains. We find that, even though the reaction rates in our reaction-diffusion model are constant, the occurrence probabilities of the receptor-scaffold reactions considered here can be strongly inhomogeneous across synaptic domains. In particular, we find that the occurrence probabilities of reactions decreasing/increasing the scaffold number, and increasing the receptor number, trace the approximate domain profile, with maxima close to the domain center. In contrast the occurrence probability of reactions decreasing the receptor number is minimal at the domain center, and shows local maxima just outside the synaptic domain. Thus, our model predicts that receptors tend to be removed from the cell membrane in membrane regions adjacent to synaptic domains, which is consistent with experimental observations Blanpied2002 .

As discussed in Sec. V, our stochastic lattice model provides a simple explanation for spatial heterogeneity in the occurrence probabilities of receptor and scaffold reactions across synaptic domains Czondor2012 ; Blanpied2002 ; Earnshaw2006 in terms of steric constraints on the receptor and scaffold membrane patch occupancies, together with the reaction-diffusion instability Haselwandter2011 ; Haselwandter2015 ; Turing1952 of the model discussed here. Similarly, we have shown previously Kahraman2016 that our stochastic lattice model provides a simple physical mechanism for distinct receptor and scaffold lifetimes at the membrane inside and outside synaptic domains Czondor2012 ; Blanpied2002 ; Earnshaw2006 . We also find Kahraman2016 that, based on the reaction and diffusion properties of synaptic receptors and scaffolds suggested by previous experiments and mean-field calculations Kirsch1995 ; Meier2000 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Ehrensperger2007 ; Calamai2009 ; Haselwandter2011 ; Haselwandter2015 , the stochastic lattice model discussed here can yield the turnover times of receptor and scaffold populations observed at synaptic domains Choquet2003 ; Specht2008 ; Calamai2009 , and predicts Kahraman2016 single-molecule trajectories consistent with experimental observations Choquet2003 ; Choquet2013 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Specht2008 ; Triller2005 ; Triller2008 .

Many essential cellular processes rely on the organization of membrane proteins into membrane protein domains Lang2010 ; Simons2010 ; Rao2014 ; Recouvreux2016 . Membrane protein domains are characterized Lang2010 ; Simons2010 ; Rao2014 ; Recouvreux2016 ; Choquet2003 ; Specht2008 ; Triller2005 ; Triller2008 ; Ribrault2011 ; Choquet2013 ; Ziv2014 ; Salvatico2015 by low protein copy numbers ($\approx 10$–1000) and protein crowding. Using synaptic membrane protein domains Ziv2014 ; Salvatico2015 as a model system, we have studied here in detail a stochastic lattice model Kahraman2016 of protein reaction-diffusion processes in crowded cell membranes. Our stochastic lattice model links the molecular noise inherent in reaction-diffusion processes to collective fluctuations in synaptic domains, and allows prediction of the stochastic dynamics of individual synaptic receptors and scaffolds. We find that molecular noise can yield substantial deviations from mean-field results, and that stochastic lattice models can be employed successfully to provide quantitative insights into single-molecule and collective dynamics of membrane protein domains Lang2010 ; Simons2010 ; Rao2014 ; Recouvreux2016 ; Choquet2003 . We focused here on the most straightforward scenario of a 1D system, which already captures Kahraman2016 the basic phenomenology of the observed fluctuations at synaptic domains. Generalization of our KMC simulations to 2D systems, which can be handled efficiently Kahraman2016 with the computational approach Elf2003 we use here, will allow more detailed and quantitative model predictions pertaining to, for instance, the stochastic trajectories of individual synaptic receptors and scaffolds at the cell membrane, which can now be directly measured Hell1994 ; Choquet2003 ; Huang2010 ; Kusumi2014 in experiments Choquet2003 ; Choquet2013 ; Meier2001 ; Borgdorff2002 ; Dahan2003 ; Hanus2006 ; Specht2008 ; Triller2005 ; Triller2008 . Our work sheds light on the organizational principles linking the collective properties of biologically important supramolecular structures, such as synaptic membrane protein domains, to the stochastic dynamics that rule their molecular components.

Acknowledgements.

We thank R. A. da Silveira, M. Kardar, and A. Triller for helpful discussions. This work was supported by NSF award numbers DMR-1554716 and DMR-1206332, an Alfred P. Sloan Research Fellowship in Physics, the James H. Zumberge Faculty Research and Innovation Fund at USC, and the USC Center for High-Performance Computing. We also acknowledge support through the Kavli Institute for Theoretical Physics, Santa Barbara, via NSF award number PHY-1125915.

Appendix A Probability distribution of jump times in a chain of Poisson processes

For completeness, we compute here the probability distribution of jump times between arbitrary states in the chain of Poisson processes in Eq. (35). We first note that

[TABLE]

Taking the Laplace transform

[TABLE]

the convolution in Eq. (60) can be conveniently represented in the form

[TABLE]

From Eq. (39) we have

[TABLE]

which implies that

[TABLE]

Transforming Eq. (64) back by using an inverse Laplace transform, one obtains

[TABLE]

Following similar steps as outlined above, the probability distribution of jump times for an $n$-step jump can be expressed as an $n$-fold convolution, resulting in

[TABLE]

which provides the probability distribution of jump times between arbitrary states in the chain of Poisson processes in Eq. (35).

Appendix B Mean jump time in a chain of Poisson processes

In this appendix we formally derive the expression in Eq. (40) for the average jump time from a state $p$ to a state $q$, $\langle t\rangle_{p\to q}$, in the chain of Poisson processes in Eq. (35). Consider the mean jump time from state $i$ to state $i+2$ in Eq. (35):

[TABLE]

where we used that $\int_{0}^{\infty}dtP_{j\to j+1}(t)=1$. Decomposing an arbitrarily long sequence of jump processes into pairs of jump processes and using Eq. (B), Eq. (40) follows from Eq. (39). Using Eq. (66), Eq. (40) can also be derived directly from the probability distribution of jump times for an $n$-step jump in Eq. (35).