Cooperative stochastic binding and unbinding explain synaptic size dynamics and statistics

Synapses are dynamic molecular assemblies whose sizes fluctuate significantly over time-scales of hours and days. In the current study, we examined the possibility that the spontaneous microscopic dynamics exhibited by synaptic molecules can explain the macroscopic size fluctuations of individual synapses and the statistical properties of synaptic populations. We present a mesoscopic model, which ties the two levels. Its basic premise is that synaptic size fluctuations reflect cooperative assimilation and removal of molecules at a patch of postsynaptic membrane. The introduction of cooperativity to both assimilation and removal in a stochastic biophysical model of these processes, gives rise to features qualitatively similar to those measured experimentally: nanoclusters of synaptic scaffolds, fluctuations in synaptic sizes, skewed, stable size distributions and their scaling in response to perturbations. Our model thus points to the potentially fundamental role of cooperativity in dictating synaptic remodeling dynamics and offers a conceptual understanding of these dynamics in terms of central microscopic features and processes.

1 Stability analysis for bidirectional cooperativity: mean-field approximation In order to determine the appropriate parameter regime for the bidirectional cooperative model, we first analyze its rate equations in the mean-field approximation. This approximation neglects fluctuations altogether, but can provide some information on qualitative behavior on e.g. stability and the existence of phase transitions. We define Θ = S/M to be the fraction of occupied sites on the lattice; with this notation, the fraction obeys the average equation Eq. (S1) has two stationary states, Θ * = 1 and Θ * = −α/(λ on − λ off ). While the full matrix (Θ * = 1) state is always a solution of this equation, some parameters combinations render the latter fixed point outside the admissible interval between 0 and 1. A linear stability analysis reveals that actually only one of these fixed points is stable -they switch their respective stabilities when the expression −α/(λ on − λ off ) enters the interval 0; 1 (see Fig. SI1).

A stochastic model with global cooperativity
In the full stochastic bidirectional cooperativity model (probabilities defined by Eq. 4 in the main text), binding and unbinding depends on the number of occupied nearest-neighbor lattice sites. This is biophysically plausible since inter-molecular interactions are of limited range. A mathematically simpler model can be defined where binding and unbinding are stochastic events, but their probabilities exhibit a global bidirectional cooperativity: they depend on the entire occupancy of the lattice, regardless of where the occupied sites are located: k on = λ on Θ + α , and k off = λ off (1 − Θ) .
Below we present numerical simulations of this model and analysis of its Master equation for the probability distribution of synaptic size.

Statistical properties of the global cooperativity model
Simulations of this global interaction model reveal that it displays statistical properties similar to the local bidirectional model. Fig. SI2 shows the analog of Fig. 4 in the main text, demonstrating the shape, stability and scaling properties of synaptic size distributions in the global interaction model. Fig. SI3A-D displays properties of the temporal fluctuations in synaptic size in this model, in analogy to Fig. 5 in the main text. The temporal properties are similar. However, the spatial patterns consisting of clusters shown in Fig. 6 of the main text are completely absent, as clearly expected, since interactions are homogeneous throughout the lattice (Fig. SI3E). The similarity of distribution shape and scaling properties between the global and local cooperativity model motivated us to investigate the Master Equation for synaptic size in the global model in more detail. Since all probabilities are determined by one variable -total synaptic size -the equation is relatively simple and can be analyzed to help understand the origin of skewed distributions in the presence of bidirectional cooperativity.

Master Equation for global models -general
In order to analyze the distribution of synaptic sizes, we formulate the master equation for stochastic binding and unbinding of molecules to the membrane matrix, where cooperativity may be introduced globally. We formulate this as a general model with different parameter sets corresponding to the Langmuir, Contact or global bidirectional cooperativity model. All these models share the property that binding and unbinding can (at most) depend on the total binding state of the synapse and not on its spatial configuration. In the following, we introduce the abbreviations A(Θ) = k on (1 − Θ) and B(Θ) = k off Θ for the total binding and unbinding fluxes per unit time, respectively. This notation is general for all models discussed in the main text. These models differ in their dependence of binding/unbinding probabilities on the fraction of bound synaptic molecules; specifically, Langmuir model k on = α , k off = β , Contact process k on = λ on Θ + α , k off = β , Bidirectional cooperative k on = λ on Θ + α , k off = λ off (1 − Θ) .

(S3)
In the most general case a non-cooperative unbinding rate β can be added also to the bidirectional cooperativity model. These in turn define the specific functions A(Θ) and B(Θ) for each model.
In a discrete master equation for a population of synapses, one only considers changes of the system size S by binding or unbinding of a single molecule; the time evolution of the probability P[Θ, t] is then In principle, the probabilities P[Θ, t] are only defined on a the discrete set of values of Θ with distance 1/M to each other. For system size M large enough, however, we can assume Θ to be (almost) continuous. This allows to approximate the discrete coupling to neighboring probabilities in Eq. (S4) via series expansions in orders of 1/M . Zeroth-order terms of this expansion, A(Θ)P[Θ, t] and B(Θ)P[Θ, t], cancel with their counter-terms in the original master equation (S4). After keeping only terms up to second order in 1/M , we arrive at a Fokker-Planck equation, as an approximation to the original master equation, Eq. (S4). The stationary distribution, ∂ t P [Θ] = 0, to Eq. (S5) is given by the expression where the constant c ensures normalization of the stationary distribution, Similar to the introduction of the non-cooperative binding rate α, that prevents the dynamical system collapsing to its empty state, a non-cooperative unbinding rate β is needed to prevent a singularity at the full state. For the simulations presented in the main text this addition is not necessary -any practical simulation time will never lead even near such a full state, if not started in its vicinity. Measured fluctuations in simulations cannot sample this tail of the distribution, as we will see in the ensuing results: distributions of synaptic size decay fast enough in their variable Θ. However, the stationary distributions in Eq. (S6) effectively includes an infinite time limit and thus requires at least an infinitesimally small β. This can also be observed by inspecting the definition of the unbinding rate, where in the latter expression the difference for small β is negligible. We are aiming to explore the effects of cooperative rates, which need to be much larger than non-cooperative rates in order to lead to the skewed distributions we want to explain.
Stationary synaptic size distributions can be calculated for each one of the three models using (S6). We compared the solutions arising from the Fokker-Planck equation with those obtained from stochastic simulations, for several cases of the different models, as shown in Fig. SI4. Both methods show a Gaussian-like symmetric distribution for the Langmuir model and for the contact process in the super-critical regime. Both exhibit a slightly skewed distribution with mean near zero for the contact process in the sub-critical regime with weak spontaneous binding. Finally, both show also similar results for the bidirectional cooperativity model. This agreement serves as a control that validates the results of the simulations shown the main text. In particular, it indicates that the skewness obtained for the bidirectional cooperativity model stems from the model itself rather than from irrelevant factors of the simulations such as boundary effects.
To analyze the stationary distribution P in more detail, we use the most general probabilities per unit time k on = λ on Θ + α, k off = λ off (1 − Θ) + β, allowing to treat all models simultaneously. The various models emerge as different limits of parameter values of this general form (compare to Eqs. (S3)). To proceed analytically, we carry out the integral in the first term of (S6). Its general structure contains three terms, where Q(Θ) and K(Θ) are low-order polynomials in Θ. These polynomials are specified for the different models in Table 1. Depending if one (or more) of the four parameters (λ on and λ off , α, or β) are zero, some of these terms in (S7) may vanish. These expressions allow us to identify relevant parameters combinations, which will determine the moments and shape of the steady state distributions. These are the averages and differences of the cooperative and non-cooperative binding and unbinding, respectively:  Table 1: Coefficients and polynomials for the solution of the integral in (S7). Only parameters with a values of zero are indicated in the first four columns, without an entry these parameters are assumed to be larger than zero.
These are the combinations which appear in the analytic solutions of the steady-state distribution. Table  1 below lists the solutions to the integral (S7) for all possible combinations of neglecting zero, one, or two of the original parameters, already using this transformation. These combinations are therefore natural candidates to be used to characterize the different behaviors in parameter space.

Analysis of the bidirectional cooperativity model using the Fokker-Planck solution
Moments in parameter space. Using the closed-form solution formulated in Eqs. (S6), (S7), and using the specific form for the bidirectional cooperativity model, one may scan the space of parameters efficiently and compute the moments of the steady-state distribution. Fig. SI5 shows the variance and skewness of the steady-state solution in a two-dimensional projection in parameter space. We choose two dimensionless parameter combinations based on the general solutions in table 1 and on the results of stochastic simulations presented in the main text. One direction for the projection (vertical axis in the figures) is µ 2 /µ 1 = (α + β)/(λ on + λ off ). This parameter combination characterizes the general strength of non-cooperative binding relative to cooperative binding. This parameters affects the skewness significantly; a decrease in skewness is seen as a function of this parameter (going up in the vertical direction for all values on the horizontal axis). This is in agreement with our general conclusion on the importance of cooperative binding in generating skewed distributions; given the reasonable assumption that both cooperative and some degree of non-cooperative binding and unbinding exist, this highlights the necessity of dominance of the first over the second for skewed distributions to appear.
The second direction is −δ 1 /µ 1 = (λ off − λ on )/(λ off + λ on ); this dimensionless parameter represents the difference between cooperative binding parameters, normalized by the total strength of cooperative binding. In agreement with the more restricted parameter range spanned by the stochastic simulations, the calculation reveals that the difference between the cooperative parameters affects the skewness slightly and only very close to the phase transition a sharp decrease in skewness is observed (Fig. SI5B).

2). (A)
For a small cooperativity ratio µ 1 /µ 2 (non-cooperative rates are much larger than cooperative rates), we find an almost constant variance of the stationary synaptic size distribution, and its shape is close to Gaussian. This fact can also be seen in the very small values of skewness shown in (B). If cooperative rates are large, the actual value depends non-monotonically on various parameter combinations. In general, skewed distributions are obtained when cooperative are larger than non-cooperative rates (see panel (B)). Both panels use δ 2 = 0, such that both non-cooperative rates are assumed identical, α = β. Relaxing this assumption, we find more complicated (non-monotonic) dependence for the exact value for the moments, although the general trend is similar to the shown behavior. The green line indicates equality in the condition given by Eq. (S20), skewed distributions are found below this line. The black line shows the parameter scan presented in Fig. 7 of the main text -deviations from the roughly Gaussian distribution are clearly visible below the crossing of these two colored lines.
Analytic condition for skewed distributions. For the full set of parameters (λ on , λ off , α, β > 0), the polynomial Q(Θ) in (S7) is given by Q = A + B. Thus, we can write the stationary probability distribution of occupied sites as product of three terms, which arise from the three terms in the integral computed above, (S7), while c is the normalization constant (c = exp(c ) from Eq. (S6) above). Here, the expression for BS(Θ) also contains the logarithmic term in Eq. (S6). Explicitly, we have In the following, we analyze the shape and characteristic length scales for each of those terms. It will be seen that BS(Θ) ("Bell-Shaped") is a symmetric Gaussian-like curve, whereas E(Θ) is exponential-like and AT (Θ) a skewed function (originating from the arctanh expression). Therefore, within the product of symmetric and skewed function, relative length-scales are important, and in particular the shape of the product will be dominated by the factor with the shortest length-scale. In this context, we define length scales as a change in Θ, which also changes the value of these functions significantly. This allows to check if any of these terms is dominant and determines the shape of the stationary distribution P [Θ] on its own. Each term is treated separately in a paragraph below. First, we inspect the central term BS(Θ) in more detail. We find that A + B is a parabola curving downward with a single maximum. Moreover, for any realistic parameter combination, we find A + B > 0 for all Θ ∈ 0; 1 , such that (A + B) can be "normalized" and seen as a distribution itself. Including the exponent 2M C log − 1, we find BS(Θ) is symmetric with respect to its mode, and looks in general bell-shaped. The mode Θ md of this bell-shaped curve (A + B) 2M C log −1 is given by Its value is important for estimating scales of the other terms in the full distribution P [Θ]: we are interested in how the slopes of this bell-shape are changed near its peak. Knowing the general shape of BS(Θ), we estimate its length scale by using the variance of the peak, where m = 2M C log − 1. For m < 0, the bell-shaped parabola flips to an U-shaped curve, shifting almost all mass of (A + B) m to the boundaries. In this case, the "variance" is constant and close to 1 as the U-shape spans the entire interval -we find that (S11) is independent of m. More relevant is the case m > 0: Using m as an auxiliary variable, find that Var (A + B) m ≈ Var A + B /m holds to a very good extent (up to a numerical constant of ≈ 0.8). For large system sizes M (above a certain critical size M crit ) we know that m and M are linearly dependent (see Fig. SI6) and thus Var (A + B) m ≈ Var A + B / M/M crit . In particular, at m = 1 the scaling between the variances is the trivial identity. Thus, we can set 1 = m using m = (2M crit C log − 1), from which we find M crit = C −1 log . The length scale L BS of the bell-shaped part BS(Θ) is hence In addition, we can carry out the integrals dΘ Θ k (A + B) with k = 0, 1, 2 to obtain the variance of the "distribution" The nonlinear term AT (Θ) generated by exponentiating the arctanh-expression in Eq. (S7) is very close to exponential at the mode of the bell-shaped term, but significantly deviates at the boundaries.
Combining these expressions, we find that the length scale (ie. width) of the symmetric, bell-shaped part BS(Θ) is given by (S13) Next we check the second term in the product E(Θ), given by a simple exponential, cf. (S9a). Thus, the length scale L E can be directly read from the expression itself, which is of the form E(Θ) ∼ exp Θ/L E . Consequently, we arrive at L E = 1/ 2M C lin , or, using the explicit expression for C lin in Tbl. 1, For the last term in the product, AT (Θ), numerical evaluation shows that it can change its value quite significantly over the range of Θ. However, it is approximately exponential for values of Θ close to Θ md given by Eq. (S10) (see Fig. SI6C). Thus, we generalize the length scale definition from the previous paragraph on E(Θ) by setting which again assumes that in the vicinity of Θ md we have AT (Θ) ∼ exp Θ/L AT . Notice that this expression generates exactly the length scale L E in (S14) when applied to E(Θ) in (S9a). After inserting the expression for AT (Θ) from (S9c) we obtain Referencing the results listed in Tbl. 1 and inserting the mode Θ md = 1 2 − δ2 2µ1 , cf. (S10), we get 2µ 1 δ 2 2 + 2µ 3 1 + µ 2 1 µ 2 δ 1 δ 2 2 + δ 1 µ 1 µ 2 − δ 2 µ 2 1 − δ 2 µ 1 µ 2 . (S16) Summarizing the length-scale analysis, we found that BS(Θ) is a symmetric, bell-shaped curve centered around Θ md . Moreover, both E(Θ) and AT (Θ) introduce an exponential skew to the stationary distribution of synaptic sizes. Thus, if the latter terms dominate, we will find highly skewed distributions, while as long as BS(Θ) dominates, we see an almost Gaussian distribution of synaptic sizes.
The dominant term in the stationary solution corresponds to the smallest length scale. However, in order to properly compare symmetric and skewed parts of the distribution, we have to combine the two terms E(Θ) and AT (Θ). To this end, note that applying the length scale definition of (S15) to the product E(Θ)AT (Θ) leads to the relation L −1 E,AT = L −1 E + L −1 AT . Evaluating this relation yields (S17) Now we can compare L BS and L E,AT . The smaller of these scales will dominate the stationary solution P [Θ]. If L BS is smaller we have a roughly Gaussian shape. Contrary, if L E,AT is smaller, we find P [Θ] to be a very fast decaying distribution, with a large skew. The direct comparison of the full expressions is rather cumbersome to analyze, L E,AT L BS ⇔ (S18) 1 2M µ 1 µ 2 1 + µ 1 µ 2 + δ 2 2 δ 1 2µ 2 1 + 5µ 1 µ 2 + δ 2 5δ 1 δ 2 − 4µ 2 1 − 4µ 1 µ 2 1 4 √ M µ 2 1 (µ 1 + 3µ 2 )(µ 1 + 5µ 2 ) − 5δ 2 2 µ 1 µ 2 − δ 1 δ 2 µ 1 + 3µ 2 2 In order to gain more insight, we chose a particularly simple limit of (S18). We set δ 2 = 0, such that the two non-cooperative rates α and β are equal (or at least very similar and their distance is small). The two length scales become (ignoring O 1 numerical constants): A few conclusions can be drawn by comparing (S19). When the system is large, M 1, we will usually arrive in a regime where the distribution of synaptic sizes is skewed. If the cooperative rates, λ on and λ off , are a lot larger than the non-cooperative rates, α and β, we also expect rather skewed distributions. Conversely, very similar cooperative rate constants, λ on and λ off , favor Gaussian shaped distribution.
With the scalings in Eq. (S19) at hand, we can also reevaluate the condition for skewed distributions, Eq. (S18). In particular, the reduced inequality becomes |λ on − λ off | λ on + λ off When this inequality does not hold, we expect Gaussian distribution of synaptic sizes. With this condition, Eq. (S20), at hand we can revisit results derived earlier, in particular compare it to the condition obtained from the mean field analysis: There, we derived α ≤ λ off − λ off (without including β yet) as indicating parameter values, where the dynamics will does not run into its absorbing state of a full scaffold. Using the default value M ∼ O(10 3 ) for the system size, which is roughly equal to µ −1 2 = 1/(α + β), both of these conditions actually restrict parameters for stable and skewed distributions to similar ranges of parameters space: Adding large enough cooperative rates is already sufficient to explain observed synaptic size dynamics.