SUPPORTING INFORMATION Mutual Repression enhances ... - PLOS
1
SUPPORTING INFORMATION Mutual Repression enhances the Steepness and Precision of Gene Expression Boundaries Thomas R. Sokolowski1 , Thorsten Erdmann2 , Pieter Rein ten Wolde3,∗ 1 FOM Institute AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands 2 University of Heidelberg, Institute for Theoretical Physics, Philosophenweg 19, 69120 Heidelberg, Germany 3 FOM Institute AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands ∗ E-mail: [email protected]
1 1.1
Methodic details Measurement of the boundary width
By default we determine the boundary width in the following two ways: Let cm n,s be the copy number of species s in a nucleus with angular index m < Nφ and axial index n < Nx , where Nφ is the number of rows around the circumference of the cylinder, and Nx is the number of colums in the axial direction along the AP axis. To compute the boundary width of the expression m m = (cm domain of a gap protein s, we compute for each row m Tn,s n,s − θs ) · (cn+1,s − θs ) as a function of n, 1 m where θs is half the copy number expected at full activation. A boundary position xm t = x (nt + 2 ) is m m defined as the position (nucleus) where Tnt ,s < 0. The values of xt are recorded in a histogram; here, the positions for the different rows m are put in the same histogram. The histogram is normalized at the end of the simulation, and the boundary width ∆x is calculated as the standard deviation of this histogram. 0 Secondly, at the end of the simulation, the slope of the average, hH(xt )i , and the standard deviation of the total Hb copy number σH (xt ) at the Hb boundary position xt are calculated from the timeσ(xt ) and φ-averaged profiles. From this, an approximation for the boundary width given by ∆x ≈ |hH(x 0 t )i | is obtained, following [1–3]. To this end, first xt is determined in the same way as in the runtime measurements, only now working on the (both time- and circumference-) averaged profile. We describe 0 in the following section how the steepness hH(xt )i is measured.
1.2
Measurement of the profile steepness
In our discrete system the measurement of a local derivative at the boundary position xt is a process prone to even small stochastic variations if a naive measurement technique is chosen. If the average boundary position xt for a set of different samples with identical initial conditions always is in between two particular nuclear positions x(n0 ) and x(n0 + 1), then using linear differences to determine the 0 steepness hH(xt )i at the boundary position may give a reasonable estimate. If, however, xt fluctuates around a particular nuclear position x(n0 ) among different samples and hH(x(n0 − 1))i − hH(x(n0 ))i significantly differs from hH(x(n0 ))i − hH(x(n0 + 1))i, the linear differences method will produce a large error bar and also markedly affect the mean of xt among these samples. As a result both the measured steepness and the quality of that measurement for a given set of parameters depends on whether xt accidently happens to predominantly vary in the interval between the same nuclear positions or not. To overcome this illness we measure the boundary steepness from the average protein profile by a two-step polynomial fitting procedure: First we fit a polynomial of 3rd degree to a region of the data around xt that contains at least four points (nuclei). The derivative of the polynomial at xt gives an initial estimate of the boundary slope, which we use this to calculate the approximative x-interval over which the profile falls from maximal to minimial expression level. If the latter is larger than the original fitting range (which usually is the case) we repeat the fitting on the enlarged interval. Since the profiles to a good
2
approximation are sigmoidal functions this improves the quality of the fit. The measured boundary slope then is defined as the derivative of the polynomial function at xt after the second fitting.
1.3
Number of cortical nuclei at cell cycle 14
The development of the Drosophila embryonic syncytium starts with a single nucleus. The first 9 nuclear divisions happen in the yolk. During cell cycles 7 to 10 a migration of the nuclei towards the cortex can be observed. However, approximately 200 polyploid nuclei stay behind in the yolk and stop dividing after their 10th cycle [4]. This quiescence persists during subsequent cell cycles, including cycle 14. As an effect of this, the number of nuclei at the cortex in cycle 14 is considerably lower than 213 = 8192. An estimate of the reduced number of cortical nuclei is given by: Ncortex ' (29 − 200) · 24 = 4992
(1)
This number indeed is closer to 212 than to 213 . Note that in our model the precise number of nuclei does not matter, rather it is the distance between the nuclear compartments and the diffusion correlation length that impact on the results. Our values for both the internuclear distance and the nuclear diameter correspond to the experimental values reported by Gregor et al. [2, 5].
1.4
Predicted copy numbers and effective protein lifetime
Our main observables are the total copy numbers of Hb and Kni, defined as follows: m m H ≡ cm n,H = cn,HM + 2cn,HD + 2
5 X
cm n,K1 j
j=0 m m K ≡ cm n,K = cn,KM + 2cn,KD + 2
5 X
cm n,H1
(2)
j
j=0
Here, for G ∈ {H, K}, cm = 1 if the promoter of species G is binding j morphogen molecules and one n,G1 j
can be equal to one for only one (repressing) gap dimer; evidently, at any given moment in time cm n,G1 j
j ∈ {0..5}. The ratio between the number of monomeric and the number of dimeric proteins is a nontrivial function of the monomer production rate, the monomer and dimer degradation rates and the parameters that determine the dimerization and dedimerization reactions. To obtain an estimate for the expected copy numbers of monomers and dimers of gene g we solved the mean-field rate equations for a simplified model which comprises monomer production, (de)dimerization and monomer and dimer degradation only, i.e. in which promoter state fluctuations and diffusion are neglected, in the steady state. We assume here that stochastic monomer production events can be accounted for by an effective mean-field production rate hβi = βhH50 i for Hb and similarly for Kni, which depends on promoter (un)binding parameters and the particular morphogen and repressor levels. This yields the following prediction for the copy number
3
of monomers (GM,hβi ) and dimers (GD,hβi ): ( 1 D D 2kon µD − koff µM − µM µD GM,hβi = D 4kon µD ) q 2 D D Dµ D + 8hβikon D koff + µD + µM koff + µD − 2kon µD GD,hβi =
1 D µ2 8kon D
( D 2 D D koff µM + µD 4hβikon + µM µM − 2kon − µM
q
) D Dµ 8hβikon D koff + µD
2 D +µ D + µM koff D − 2kon µD
(3)
D D (koff ) are the dimerization forward Here µM (µD ) is the monomeric (dimeric) degradation rate and kon (backward) rates, respectively. From this we calculate the total expected copy number Ghβi := 2GD,hβi + GM,hβi at effective production rate hβi. In particular in the full-activation case, i.e. when the probability to be fully activated and unrepressed hG05 i ≈ 1 and therefore hβi ≈ β, the above estimates correspond to average values from our simulations very well. 1 (µM GM + 2µD GD ). Our stanWe define the effective degradation as µeff = µeff (GM , GD ) = GM +2G D −3 dard values result in µeff ≈ 4.34 · 10 /s with GM = GM,β and GD = GD,β .
2 2.1
Additional analysis Poissonian limit with dimerization
In [3], it was shown that, when D → ∞, the variance in the protein concentration becomes equal to the mean concentration: diffusion washes out bursts in gene expression, thus reducing the non-Poissonian part of the noise. However, in that model the proteins do not dimerize, in contrast to our model. With dimerization, a different limit for the variance in the total protein concentration is approached as D → ∞. To derive this limit, first note that the total protein copy number G of a protein G is G ≈ 2GD + GM . Assuming that hGD GM i ≈ hGD ihGM i (our simulations indicate that this approximation is very accurate), 2 we find that the variance σG in G is: 2 2 2 σG ≈ 4σG + σG , D M
(4)
2 2 where σG is the variance in the dimer level GD and σG is the variance in the monomer level GM . Both D M monomers and dimers are subject to spatial averaging, and therefore their variances can be written in the form [3]:
1 2 σ0,G − GM M N 1 2 = GD + σ0,G − GD D N
2 σG = GM + M 2 σG D
(5)
2 Here N is the number of nuclei contributing to the averaging, which is proportional to D, and σ0,G is M/D the variance in the monomer and dimer levels in the absence of diffusion, respectively. The part preceded by 1/N represents the variance that can be reduced by spatial averaging. Plugging these expressions into
4
the previous and using G = 2GD + GM we arrive at: 1 2 2 2 4σ0,GD + σ0,G − 4GD − GM σG = 4GD + GM + M N 2GD 1 2GD 2 = 1+ G+ σ0,G − 1+ G G N G 1 2 σ =: (1 + fD ) G + − (1 + fD ) G N 0,G
(6)
Note that N is the same for both monomers and dimers because their diffusion√constant does not differ in our p model. Evidently, the lower bound for σG in the limit N → ∞ is not G any more, but given by (1 + fD ) G, where fD is the fraction of proteins in the dimer state with respect to the total protein number (implying fD ≤ 1). This is indeed what we observe in our data for σG . In our simulations the equilibrium is strongly shifted towards the dimerized state, so that fD ≈ 0.97. We can understand the limit N, D → ∞ intuitively by noting that in this limit there is no noise in the nuclear protein concentration due to the stochastic production and decay of molecules in each of the nuclei—this is because the synthesized molecules are immediately donated to a reservoir that is infinitely large; instead, there is only noise in the nuclear protein concentration due to the sampling of molecules from this 2 2 reservoir, which obeys Poissonian statistics: σG = GM and σG = GD . This yields, for N, D → ∞, M D 2 2 2 σG = 4σGD + σGM = 4GD + GM = (1 + fD )G.
2.2
Bifurcation analysis
In order to predict the regions in which bistability can be expected for different amplitudes A of the morphogen gradients we performed a deterministic mean-field bifurcation analysis for a simplified 1dimensional version of our model of mutual repression between hb and kni . The analysis is based on the following two equations describing the change of the mean-field total copy number of Hb (H(x)) and Kni (K(x)) at position x: ∂t H(x) = βH (x) ∂t K(x) = βK (x)
2 KR 2
− µH H(x)
(7)
2
− µK K(x)
(8)
2 + [f K(x)] KR D 2 KR
2 + [f H(x)] KR D
Here βH and βK represent the protein synthesis rates, µH and µK the corresponding (effective) degradation rates, KR is the dissociation constant of cooperative repressor binding to the promoter and fD is the fraction of proteins in the dimerized qstate. Note that, since the intermediate step of dimerization is
R /k R if k R and k R are the binding rates of the dimers. To neglected here, we have to take KR = koff on on off facilitate calculations we make two further simplifying assumptions here:
1. We neglect activation dynamics and resulting promoter state fluctuations, i.e. we assume that certain constant levels of the activators at position x lead to average constant production rates βH = β([Bcd](x)) and βK = β([Cad](x)), respectively. In our standard case β([Act](x)) = [Act]5 (x)/([Act]5 (x) + 6905 ) for both [Act] = [Bcd] and [Act] = [Cad]. 2. In our simulations we have different degradation rates for monomers and dimers so that the effective total degradation rate depends on the monomer-to-dimer ratio, which in turn varies with the total copy number (see section 1.4). Thus, in principle, also fD and µH and µK are functions of x, or the corresponding activator levels. Since this introduces further nonlinearities into the above equations and complicates their solution, we substitute the degradation rates µH and µK by a
5
constant value µeff , which is the effective degradation rate for the maximal expression level (full activation). Also for fD we take the constant value for full activation, fD ' 0.97, which reflects that the dimerization equilibrium in our simulations is strongly shifted towards the dimerized state. The predictions concerning the bifurcation behavior only change marginally if µeff and fD values for lower expression levels are used. For each position x with local activator levels corresponding to the ones in the simulations we calculated fixed point solutions for the copy number pair (H(x), K(x)) starting from the steady-state assumption ∂t (H(x), K(x)) = (0, 0). The stability of the fixed points was determined starting from the Jacobian for the above ODE system: ∂H [∂t H] ∂K [∂t H] J(H, K) = (9) ∂H [∂t K] ∂K [∂t K] Within the relevant parameter regime we obtained fixed points with either two negative eigenvalues (i.e. stable fixed points) or one positive and one negative eigenvalue (i.e. saddle points). The determinant therefore completely characterizes the stability of the fixed points. If det J(H0 , K0 ) < 0, then (H0 , K0 ) is a saddle point. Otherwise it is stable. Fig. S1 shows the fixed point solutions for Hb and Kni as a function of x for different activator amplitudes A. Stable solutions are drawn with solid, unstable solutions with dashed lines. Depending on the A value, the system displays a saddle node bifurcation at a point towards the anterior (Hb) or posterior (Kni) from midembryo. Within the region confined by the bifurcation points two stable and one unstable fixed points exist for each gene, implying bistability. The region clearly widens for increasing A and spans almost the whole embryo length for A = 8. Our deterministic analysis therefore predicts the enlargement of the region of bistability as observed in our single nucleus simulations.
2.3
Estimation of switching times
To quantify the swiching times in the presence of bistability we performed simulations of isolated single nuclei featuring the same set of reactions and parameters as in the full scale simulation. To obtain estimates of switching times at different positions x along the AP axis we set the levels of Bcd and Cad in the given nucleus equal to the ones at x in the space-resolved simulations. The switching time was estimated by calculating from long time trajectories of the total Hb and Kni copy numbers the relaxation time ts of the average correlation function hC(t)it0 ≡
hIH (t0 )IK (t)it0 hIH (t0 )it0
(10)
where IH (IK ) are indicator functions which are one if the difference in the total gap gene copy numbers ∆N = H − K is above (below) a certain threshold ΘN (−ΘN ). ΘN thus defines the regions within which the switch is considered to have switched to the Hb–high or Kni–high states, respectively, and serves to separate the stable attractor states from the transition region. We found that ΘN = 200 is a reasonable choice for our set of parameters. We determined the switching times from one long sample for different positions x and different activator amplitudes A and find that ts is very similar within the double-activated bistable regions for high A. To obtain an error estimate we additionally calculated block averages of estimated switching times among 10 long samples for various A at midembryo (x = L/2). Table S1 shows our results from the latter procedure. Note that for A = 1 the system is not truly bistable yet because for A = 1 we have half-activation at midembryo and due to the lack of diffusion large promoter-state fluctuations dominate over long-time switching potentially induced by mutual repression. Consequently, the given number does not reflect a switching time. We cite it here for completeness, however.
6
Table S1. Switching times at midembryo for different activator levels. Activator amplitude A Switching time ts [s] (1) (6343.7 ± 17.2) 2 20302.5 ± 74.5 4 20957.3 ± 54.9 8 20994.7 ± 67.2
2.4
Analysis of statistical properties of the boundary
Our measures for both the boundary steepness and the variance of the boundary position are based on averages over both the time and the circumference of the embryo which were calculated during runtime. While the double-averaging procedure limits the amount of data that must be stored and facilitates rapid acquisition of good statistics, it also discards information about the microscopic properties of the boundary at a given time instance. Based on the average data it is impossible to determine whether the blurring of the boundary quantified by ∆x is due to concerted stochastic movements of a steep and rather homogeneous instantaneous boundary or simply due to stochastic fluctuations of the boundary position in each nuclear row around a well-defined constant mean boundary position (or due to both). In the latter case the boundary will be rough at each given time instance, i.e. the time average of the boundary position variance in the cicumferential direction will be large, but the time variance of its circumferential mean will be negligible. The opposite will be the case in the other extreme. These quantities therefore can be used to distinguish the two hypothetical situations. The overall boundary width in both cases is given by the sum: 2 ∆x2 = σx2t (φ) + σhx t iφ (t)
(11)
Here h. . . iφ denotes the average over the circumference, while the bar denotes the time average. An identical variance decomposition can be made for the fluctuations of the Hb copy number at any position x along the AP axis. Similarly, comparing the average of the profile steepness for a particular nuclear row and time instance to the steepness of the time- and circumference average of the copy number reveals whether the steepness of the average profile is due to concerted movements of similarly shallow instantaneous profiles or due to unconcerted fluctuations of steep instantaneous profiles. In order to determine which of the portrayed blurring mechanisms is dominant in our system we performed the described variance decomposition for a set of 100 instantaneous outputs of the fully resolved 2D system in steady state, i.e. for 6400 different total Hb copy number profiles along the AP axis, for both the variances at the boundary and for the steepness at the boundary and for both the system with and without mutual repression. We focused on our standard parameter set (see Table S2) and a range of gap protein diffusion constants D. 2.4.1
At each time instance the boundary is rather rough
Fig. S2 shows for the systems with (Fig. S2A and C) and without (Fig. S2B and D) mutual repression the variance decomposition for the variance of the Hb copy number at the boundary (Fig. S2A and B) and for the variance of the boundary position xt (Fig. S2C and D) as a function of the Hb protein diffusion constant D. As a control we compare the total variances calculated from the instantaneous profiles to the variances accumulated during runtime and, in case of ∆x, to the value obtained from the 0 approximation ∆x = σH (xt )/|hH(xt )i | (note that here h. . . i is the average over both time and φ). We see a good agreement between these quantities. The plots reveal that both for σH (xt ) and ∆x the variance over the circumference at a fixed time is by far the dominant contribution to the overall variance. This implies that in our system the boundary is indeed very rough at each time point and that concerted boundary movements do not occur.
7
2.4.2
At each time instance the profiles are slightly steeper than their average
The calculation of the variance decomposition is less straightforward for the slope. In particular for low D, when spatial averaging is still inefficient, the instantaneous profiles are very ragged and the boundary threshold value typically is crossed at multiple positions along the AP axis. This makes it impossible to uniquely define an instantaneous boundary position as required to calculate the instantaneous boundary slope. In order to perform the analysis at least on a subset of the data we introduced a protocol which only takes into account instantaneous profiles with a single boundary crossing, rejecting all other profiles. For low D, however, the rejection rates rise above 90%. We therefore decided to smoothen the profiles by computing running averages between a fixed number ν of nuclei along the AP axis before the analysis. The averaging lowers the rejection rate dramatically, however it also decreases the profile steepness and therefore manipulates the observable of interest. Nevertheless we can make a qualitative statement on the base of the results obtained for only slight smoothening of the profiles (ν = 3). For simplicity and due to increased data abundance, in this analysis we used simple finite differences to determine the slope. In Fig. S3 we plot the average of the instantaneous boundary steepness for different degrees of smoothening (averaging over ν = 3, 5, 7 nuclei along the x-axis) as a function of D and compare to the steepness of the average Hb profile for the system with (S3A) and without (S3B) mutual repression. While the data for ν > 3 clearly must be considered biased by the running averages, the values for ν = 3 show that the instantaneous boundary slope on average is higher than the slope of the average profile, in particular for low diffusion constants. The variance decomposition for the boundary position xt shows that the variance of the circumference mean of the boundary position in time is very small. This implies that the steepness of the circumferenceaveraged profiles should be approximately equal to the steepness of the time- and circumference-averaged profile. As a control we therefore repeated the above analysis on the 100 φ-averaged instantaneous profiles of the same dataset. The averaging along the circumference significantly reduces the number of profiles with ambiguous boundary positions. We therefore were able to obtain reasonable estimates of the observable without pre-smoothening of the profiles (ν = 1). The results are shown in Fig. S3C for the system with mutual repression and Fig. S3D for the system without mutual repression. In the system with mutual repression the average slope of the φ-averaged profiles for ν = 1 agrees well with the slope of the both time- and φ-averaged Hb profile. In the system without mutual repression the φ-averaged profiles are slightly steeper than the average.
3
Additional simulations
3.1
Influence of the Hill coefficient
To address the influence of changing activator cooperativity on our results we performed simulations with reduced number of activator binding sites nmax . While in our model this is achieved by simply reducing the number of intermediate states between the empty promoter state and the producing promoter state, the binding parameters have to be rescaled with care to preserve the activation equilibrium at midembryo. Since we assume the activator binding rates to be diffusion limited, the necessary changes affect the A unbinding rates koff,n = a/bn . However, even when preserving the equilibrium, the freedom in the choice of these parameters allows for altering the time scale of transitions between the different activation levels. In order to rescale the rates in a unique fashion upon lowering nmax we imposed the following constraints: A 1. For all nmax the effective activator dissociation constant at midembryo KD = A1/2 = 690 is preserved, which implies that for all nmax the average activation probability at midembryo is 1/2.
2. The waiting-time distribution for the unbinding from the producing state is the same for all nmax A and, for comparison, equal to the one for the default cooperativity nmax = 5, i.e. ∀n : koff,n = max A const = koff,5 .
8
3. The off-rate reduction per subsequent activator binding is always 1/b, i.e. ∀n : b(n) = b. A A A Note that for nmax = 1 the first two conditions can be met together only if KD kon = koff,5 , which is not the case for our parameter set. We therefore restricted ourselves to nmax ∈ {2, 3, 4, 5}. For each nmax , the above constraints were used to uniquely determine the parameters a and b from the exact analytical solution for the average occupancy of the producing state, which was obtained from steadystate mean-field solutions of the chemical mass-action ODEs. Interestingly this results in only minor differences in a among the different nmax values, while the reduction per binding step 1/bnmax becomes significantly larger for lower nmax . This fact has an important implication for the noise charactistics of the different promoters: If anmax ' a = const for all nmax then the unbinding rate from the state binding (nmax − 1) activator proteins (the “highest” non-producing activator state) is given by: A A max −1) koff,(n ' a/bn(nmax = bnmax koff,5 max −1) A Since bnmax markedly increases with decreasing Hill coefficient the off-rate koff,(n for low nmax will max −1) be higher than the corresponding rate for high nmax . This will favor rapid returns to the producing state with nmax bound activator molecules for high Hill coefficients, whereas for low Hill coefficients the promoter is more likely to descent into the regime with less activator molecules bound. The fact that this is less likely for higher Hill coefficients is compensated by the fact that also the time to return to the producing state from the states binding low numbers of activator molecules on average is longer for higher nmax . Note that the mean off-time–just as the mean on-time–is the same for all nmax . In short, for the promoters with higher Hill coefficients we expect an off-time distribution with high probability weight on short off-times and a long low-probability tail for long off-times, while the distribution for lower Hill coefficients should resemble an exponential. In order to illustrate this effect we recorded long time-trajectories of the occupancy of the producing state in a single isolated nucleus close to midembryo for different nmax and without mutual repression nor diffusion. All other parameters were kept at the standard values. From these trajectories we determined the on- and off-times of the promoter and binned them into a histogram. The results are shown in Fig. S5. It can be seen that while for nmax = 2 the two distributions are exponential with approximately equal mean, the off-times distribution increasingly deviates from an exponential distribution as nmax is increased; more probability is shifted to very short off-times and very long off-times, causing the emergence of a long tail in the distribution.
3.1.1
Also for lower Hill coefficients mutual repression steepens profiles without corrupting boundary precision
The broadening of the off-times distribution is expected to result in higher output noise for high nmax as compared to low nmax . This is confirmed by the simulations of the full-scale spatially resolved system for different nmax . Fig. S6 shows σH (xt ), the average standard deviation of the total Hb copy number at 0 the boundary position xt (upper panels), the steepness |hH(xt )i | of the average Hb profile at xt (middle panels) and the boundary width ∆x (lower panels) as a function of the gap protein diffusion constant D for nmax ∈ {2, 3, 4, 5}. σH (xt ) is indeed decreasing upon lowering nmax , in particular in the regime of low diffusion constants. For higher diffusion constants the decrease is less pronounced: spatial averaging is efficient enough to lower the output noise down to the observed values irrespective of the width of the offtime distribution. The noise decreases less markedly for the systems with mutual repression. This is most likely due to the fact that lowering nmax also increases the probability of occasional repressor production beyond midembryo, which in turn increases the noise. The steepness plots reveal that, although the profiles naturally become less steep for lower nmax , the steepness in the systems with mutual repression is markedly higher than the one in the system without mutual repression. In all systems the steepness as a function of D shows a very similar behavior: Upon increasing D the steepness in the systems with mutual repression first increases towards a maximum before it rapidly decreases. Since both σH (xt ) and
9 0
|hH(xt )i | change with nmax in a similar fashion, in particular in the region around D = 1 µm2 /s, the width ∆x as a function of D also looks very similar in this region for all nmax . In all cases the profiles in the system with mutual repression are more precise and markedly steeper as compared to the system without mutual repression at a D-value which is one order of magnitude less than the optimal value in the systems without mutual repression. Therefore the basic effect observed in our simulations for nmax = 5 persists in the simulations for lower Hill coefficients. 3.1.2
Lower Hill coefficients allow for stronger morphogen level variations
Although lowering nmax in our system reduces the protein production noise it also markedly decreases the steepness of the gene activation profiles. An important implication of this is that for lower nmax the activation probability beyond midembryo increases. Lowering nmax thus is similar to increasing the activator amplitude A and, in principle, might result in the creation of a bistable region around midembryo already for lower A-values as compared to the system with nmax = 5. We analysed how the results for ∆x as a function of A for the case of correlated variations change as nmax is decreased. Fig. S7 shows ∆x(A) for nmax ∈ {2, 3, 4, 5} and D = 1.0 µm2 /s for systems with and without mutual repression. Overall, ∆x(A) is very similar for all considered nmax . For A ≤ 2 the width ∆x in the systems with mutual repression is always lower than in the systems without mutual repression. The minimal ∆x is attained at A = 1 in all cases. The main difference is in how ∆x changes with A for A > 1: The lower nmax , the slower the width increases with A. Thus, while lower Hill coefficients decrease the steepness of the profiles significantly, they may prove beneficial by extending the range over which extrinsic variations are successfully buffered without increasing intrinsic fluctuations of the boundary.
3.2
Influence of the expression level
In order to examine the influence of a changed signal-to-noise ratio on our results we performed simulations with altered production dynamics. We did this by (1) introducing bursty production, i.e. producing 10 copies of the gap protein monomer at a time with a 10 times lower production rate (β = β0 /10), and (2) by keeping the burst size at one and changing the production rate β. To preserve the binding equilibrium of the repression reaction at midembryo upon changing β we also changed the off-rate of the repressor dimers by a factor fβD , which is the ratio between the expected number of dimers at midembryo for the altered production rate β and the corresponding value for the standard production rate β0 . Note that, since in our system the copy numbers of both monomers and dimers depend on β in a nontrivial fashion (see section 1.4) the effective copy number increase typically does not correspond to the ratio β/β0 . Therefore fβD > β/β0 for β > β0 . 3.2.1
Bursty production has only a marginal influence on the boundary properties
In Fig. S8 we plot the standard deviation of the total Hb copy number at the boundary, the steepness of the total Hb copy number profile at the boundary and the boundary width ∆x as a function of D for the system with bursty production (burst size 10). There is no significant difference as compared to the system with normal production (burst size 1, compare to Fig. S6(D) or Fig. 3 of main article). For low D the production noise is marginally higher with bursty production, resulting in a slight increase of ∆x in this regime; the effect of varying D, however, is much more important. This is most likely a consequence of the fact that for the given Hill coefficient nmax = 5 promoter state fluctuations are already at a high level due to a very broad off-time distribution (see section 3.1).
10
3.2.2
Increased production rates reveal different noise scaling behavior for different regimes of the diffusion constant
Upon increasing the production rate and consequently the total copy number of the gap proteins we may expect a relative decrease in the output noise, but only if the latter is purely Poissonian. In our system √ this corresponds to the limit of high gap protein diffusion constants. In that limit, we expect σG ∝ G, where σG is the noise in the total gap gene copy number G. However, in the abscence of spatial averaging, i.e. for the limit D → 0, non-Poissonian noise prevails and the expected scaling is σG ∝ G [3]. If the copy number profile is scaled uniformly at each AP position x, which–to a good approximation–is the case in our system, we expect for the scaling of the gradient at midembryo G0 (xt ) ∝ G(xt ). The expected √ scaling for the boundary width ∆x then is ∆x ∝ 1 for low diffusion constants and ∆x ∝ 1/ G for high D. While the overall characteristics of the boundary are very similar to the system with β = β0 , a comparison roughly confirms the predicted scaling. Fig. S9 compares for Hb the standard deviation of the total copy number at the boundary (S9A), the steepness at the boundary (S9B) and the resulting boundary width (S9C) as a function of D for increased production rates β = 2β0 and β = 4β0 to the corresponding values for the sytem with production rate β/2. Thus, the values for β = 4β0 are compared to β = 2β0 and the values for β = 2β0 are compared to β0 . Blue lines mark the expected change of the quantities as predicted by the scaling relations, where the corresponding copy number increase is given by the factor f2 ≡ f2β0 = 2.22 and f4 ≡ f4β0 /f2β0 = 2.16. Here fβ ' fβD is the predicted total copy number at midembryo divided by the corresponding value for β = β0 . The plots show that while the slope ratio is roughly equal to f2 (f4 ) for all D both in the system with (green) and without (red) mutual repression, the noise ratio depends on the diffusion constant and also slightly differs for the systems with and without mutual repression. Nevertheless the predicted scaling behavior is confirmed√in both √ cases: in the low diffusion constant regime the noise ratio is roughly f2 together this leads to a boundary width ratio of one for (f4 ) and approaches f2 ( f4 ) as√D increases; √ low D which decreases towards 1/ f2 (1/ f4 ) for higher D.
3.3
Activation of both gap genes by a single gradient
In the one-morphogen gradient scenario, both hb and kni are activated by the Bcd gradient. Here, kni is activated in the same way as hb, namely by 5-step cooperative binding, but with a lower activation threshold. This results in the induction of both genes in the anterior half of the embryo up to the posterior Hb boundary and of kni in an additional region posterior to the Hb boundary. Given that hb represses kni more strongly than vice versa in the double-activated bistable region, this parameter choice will result in the formation of two neighboring domains. We chose the kni activation threshold to be lower by a factor of 1/2, which causes an offset of its half-activation point by approximately 10 nuclei (83 µm) towards R,K the posterior. We varied the protein diffusion coefficient D and koff , the off-rate of the Kni repressor dimers from the hb promoter. The rate for dissociation of the Hb dimers from the kni promoter was kept R,H at the standard value koff = 5.27 · 10−3 /s in all simulations. R,K The diffusion constant D of the gap proteins and the dissociation rate koff of Kni from the hb promoter are indeed key parameters. On the one hand, hb must repress kni more strongly than the other R,K way around, because otherwise there will be only one kni domain. On the other hand, when koff is high, then kni is only significantly expressed when D is low, because kni represses hb more weakly than vice versa, which means that low amounts of invading Hb dimers are sufficient to shut off Kni production almost completely; indeed, in this regime, kni has hardly any effect on the precision of the hb expression R,H R,K domain. We found that when koff /koff is roughly between 0.1 and 1, both hb and kni domains are √ R,H R,K formed robustly. In Fig. S4 we display the case for koff /koff = 1/ 10 ≈ 1/3. Fig. S4A shows that the maximum of the average Kni copy number is lower than that of Hb, even though for x < 60 %EL kni is essentially fully activated by Bcd (Fig. S4B). The lower maximum is due to the fact that hb represses kni more strongly than vice versa. Another point worthy of note is that the
11
fluctuations in the Kni copy number in the Kni domain are higher than those of Hb in the Hb domain (Fig. S4C). This is essentially due to the small width of the Kni domain: kni is either fully activated by Bcd yet still repressed by Hb or not repressed by Hb yet stochastically activated by Bcd. Panels D-F show, respectively, the noise in the Hb copy number at the hb expression boundary, the steepness of this boundary, and the width of this boundary, as a function of the diffusion constant D of the gap proteins. It is seen that the results are highly similar to those of the two-gradient motif. The noise in the Hb copy number at the boundary is not much affected by mutual repression (panel D), while the steepness, and consequently boundary precision, is markedly enhanced by mutual repression, especially when the diffusion constant is small. Note that while for the two-morphogen gradient scenario 0 the approximation ∆x ≈ σH (xt )/|hH(xt )i | is in very reasonable agreement with ∆x as measured from the distribution of threshold crossings p(x), here the agreement is much less. This is due to sporadic repression events in the anterior region where hb and kni are both fully activated, which leads to a long tail of p(x) extending towards the anterior pole; while p(x) in the tail is small, the fact that the tail is long does 0 markedly increase the standard deviation ∆x. Given that the approximation ∆x ≈ σH (xt )/|hH(xt )i | works so well for all the other cases, we consider this approximation, which does not suffer from sporadic but strong hb repression events in the anterior, to be more reliable. We therefore conclude that also in the one-morphogen gradient scenario, mutual repression can enhance both the steepness and the precision of gene-expression boundaries.
References 1. Tostevin F, ten Wolde PR, Howard M (2007) Fundamental limits to position determination by concentration gradients. PLoS Comput Biol 3: e78. 2. Gregor T, Tank DW, Wieschaus EF, Bialek W (2007) Probing the limits to positional information. Cell 130: 153–164. 3. Erdmann T, Howard M, ten Wolde PR (2009) Role of spatial averaging in the precision of gene expression patterns. Phys Rev Lett 103: 258101. 4. Foe VE, Alberts BM (1983) Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in drosophila embryogenesis. J Cell Sci 61: 31–70. 5. Gregor T, Wieschaus EF, McGregor AP, Bialek W, Tank DW (2007) Stability and nuclear dynamics of the bicoid morphogen gradient. Cell 130: 141–152. 6. Abu-Arish A, Porcher A, Czerwonka A, Dostatni N, Fradin C (2010) High mobility of bicoid captured by fluorescence correlation spectroscopy: Implication for the rapid establishment of its gradient. Biophys J 99: L33–L35.
Table S2. The standard set of the Parameter / Quantity Symbol Geometry Number of nuclei in axial direction x Nx Number of nuclei in circumferential direction φ Nφ Internuclear distance (lattice spacing) l - resulting embryo length L Nuclear radius R - resulting nuclear volume V Morphogen / activator gradients Standard copy number at the poles A0 - resulting copy number at half activation A 12 Morphogen diffusion length λ Activation dynamics Intranuclear activator diffusion constant DA Activator binding site target size α A - resulting activator binding rate kon A Activator unbinding rate koff,n Hill coefficient of activation nmax Gap gene dynamics Gap protein monomer production rate β0 Gap protein monomer degradation rate µM Gap protein dimer degradation rate µD D Dimerization forward rate kon D Dimerization backward rate koff - predicted total Hb (Kni) copy number at full activation and D = 0 Hβ0 (Kβ0 ) Gap protein diffusion constant D Repression dynamics R Repressing dimer binding rate kon R Repressing dimer unbinding rate koff 0.40 5.27·10−3
µm3 /s 1/s
1/s 1/s 1/s µm3 /s 1/s
3.37 3.37·10−2 3.37·10−3 0.80 5.59·10−3 775 varied
µm2 /s nm µm3 /s 1/s
µm
6720 690 119.5 3.2 10 0.40 410/6n 5
µm µm µm µm3
64 64 8.5 544 6.5 143.8
equal to 0.1 ·
A same as kon varied in some simulations − kA,5 ;
A same as 2kon
A = 4παDA diffusion limited, kon n = no. of bound activator molecules supported by [2]
supported by [6]
varied in some simulations supported by [2] supported by [2]
supported by [5] assuming spherical shape
supported by [2]
most important parameters in our model. Value Unit Remarks
12
13
Figure S1. Bifurcation analysis of a one dimensional mean-field model of mutually repressing gap genes activated by morphogen gradients. Plotted are the stable (solid lines) and unstable (dashed lines) fixed points of the copy number of Hb (colored lines) and Kni (grey lines) as a function of the AP position x as predicted by a bifurcation analysis performed on a 1D mean-field model in which hb and kni are activated cooperatively by their respective morphogens and mutually repressing each other. Different colors correspond to different fixed points. The different panels show the solutions for activator amplitudes (A) A=1, (B) A=2, (C) A=4 and (D) A=8. All other parameters values are the standard values from Table S2. Activator concentrations at position x used in the mean-field analysis correspond to the ones in the 2D stochastic simulations. Away from midembryo each gap protein level has only one stable fixed point and one of the two levels is always zero. For all A there is a region around midembryo in which the protein levels have two stable and one unstable fixed point, implying bistability. In this region the analysis predicts bistable switching between the high-Hb–low-Kni and the low-Hb–high-Kni state. For clarity we color-code the Hb fixed points only. The Kni solutions are identical to the Hb solutions mirrored with respect to midembryo.
14
Figure S2. Decomposition of variances at the boundary. (A) Decomposition of the total Hb copy number variance at the average boundary position for the system with mutual repression as a 2 function of the gap protein diffusion constant D. Plotted are: σH (xt ) the total (time- and circumference-) variance measured during runtime (RT, black), the same quantitity determined from a 2 set of 6400 instantaneous profiles (IP, red, 64 AP rows at 100 different time points), σhHi the variance φ
2 the time average of the variance of in time of the circumference average of H(xt , φ) (green) and σH H(xt , φ) over the circumference (blue). (B) The same as (A) for the system without mutual repression. (C) The same variance decomposition as in (A) for the Hb boundary position xt instead of the copy number. The black line shows the ∆x values measured as the standard deviation of the boundary position histogram accumulated during runtime (RT), the grey dashed line the corresponding values 0 determined from the approximation σH (xt )/|hH(xt )i |. (D) The same as in (C) for the system without mutual repression. In both cases, the main contribution to the total boundary variance σx2t comes from σx2t , implying that the blurring of the boundary is rather due to roughness than due to concerted boundary movements.
15
Figure S3. Microscopic properties of the boundary steepness. Panels (A) and (B) compare for different gap protein diffusion constants D the average Hb profile steepness at the boundary measured in a set of 6400 instantaneous profiles (64 AP rows at 100 different time points) to the 0 steepness of the (time- and circumference-) average of the Hb profile (|hH(xt )i |, black) for different numbers ν of neighboring data points used in calculating running averages over the instantaneous profiles for the system with (A) and without (B) mutual repression. Although for increasing ν the instantaneous profiles become less steep as a consequence of the smoothening, the values for ν = 3 indicate that the profiles at a given row and time instance are slightly steeper than the average profile. In panels (C) and (D) we show results of the same analysis performed on the 100 circumference-averages of the instantaneous profiles, again for the system with (C) and without (D) mutual repression. Here ν = 1 is the data obtained without calculating running averages (magenta). In both systems the steepness of the φ-averaged profiles agrees reasonably well with the steepness of the 0 average profile |hH(xt )i |.
16
Figure S4. The effect of mutual repression in a system where both hb and kni are activated by the Bcd gradient. See following page for description.
17
Figure S4. The effect of mutual repression in a system where both hb and kni are activated by the Bcd gradient. (A) Time- and circumference-averaged Hb (hHi, solid lines) and Kni (hKi, dashed lines) total copy-number profiles along the AP-axis for various Hb and Kni diffusion constants D. (B) AP profiles of the average standard deviation of the total gap gene copy number for Hb (σH , solid lines) and Kni (σK , dashed lines). Note that the noise in K in the Kni domain is larger than that in H in the Hb domain. (C) AP profiles of the probabilities hH50 i + hH51 i and hK50 i + hK51 i that the hb (solid lines) and kni (dashed lines) promoters have 5 Bcd molecules bound to them, respectively; in the absence of mutual repression between hb and kni , these profiles would directly determine the expression of hb and kni . (D) The noise in the Hb copy number at the hb expression boundary as a function of the Hb and Kni diffusion constant D. (E) The steepness of the hb expression boundary as a function of the diffusion constant of the gap proteins. (F) The width ∆x of the hb expression boundary as a function of the Hb and Kni diffusion constant. The grey line corresponds to 0 the approximation ∆x ≈ σH (xt )/|hH(xt )i |, which we consider to be more reliable than ∆x as measured from the distribution of threshold crossings, p(x); the latter suffers from sporadic but strong suppression events of hb by kni in the anterior, which leads to a long tail of p(x), increasing ∆x. It is seen that while mutual repression has hardly any effect on the noise in the copy number at the boundary, it does markedly enhance the steepness of the boundary, and thereby its precision. The ratio of the R,H R,K Hb–kni -promoter dissociation rate over the Kni–hb-promoter dissociation rate is koff /koff = 1/3.
18
Figure S5. On- and off-times distributions of the hb promoter for different Hill coefficients nmax in a nucleus at midembryo. The panels show normalized histograms of the times spent in the producing (n = nmax ) promoter state (“ON”, green) and of the times spent in the non-producing (n < nmax ) states (“OFF”, red) for (A) nmax = 2, (B) nmax = 3, (C) nmax = 4 and (D) nmax = 5 (standard case). It can be seen that with increasing Hill coefficient nmax the off-times distribution changes from an exponential to a non-exponential distribution with high weight on very short off-times (implying fast returns to the producing state) and a with a long tail of long off-times. Since the off-rate from the producing state is kept the same for all nmax the on-times distributions remain unaltered. The on- and off-times have been determined from long time trajectories (ttotal = 105 s) of the occupancy of the producing state with a sampling resolution of 0.5 s.
19
Figure S6. Boundary characteristics for reduced Hill coefficients nmax . See following page for description.
20
Figure S6. Boundary characteristics for reduced Hill coefficients nmax . The standard deviation of the total Hb copy number at the boundary (σH (xt ), upper panels), the gradient of the 0 average Hb total copy number gradient at the boundary (|hH(xt )i |, middle panels) and the boundary width (∆x, lower panels) as a function of the gap protein diffusion constant D for the systems with (green) and without (red) mutual repression and Hill coefficients (A) nmax = 2, (B) nmax = 3, (C) nmax = 4 and (D) nmax = 5 (standard case). Grey dashed lines are values determined from the 0 approximation ∆x = σH (xt )/|hH(xt )i |, solid lines are values calculated from the distributions of xt . Broad dashed lines are the values for D = 0. Black dotted lines mark the D-value where the boundaries are both steep and precise due to mutual repression.
21
Figure S7. The effect of changing the activator amplitude A on the boundary precision for reduced Hill coefficients nmax . Shown are the the boundary width ∆x with (green) and without (red) mutual repression as a function of ∆xA , the separation between the Hb and Kni boundaries expected in the system without mutual repression, and the corresponding activator amplitude A for Hill coefficients (A) nmax = 2, (B) nmax = 3, (C) nmax = 4 and (D) nmax = 5 (standard case). In all cases D = 1.0 µm2 /s.
22
Figure S8. The effect of bursty gap protein production on the Hb boundary precision. The plot shows σH (xt ) the standard deviation of the total Hb copy number at the boundary, the 0 steepness |hH(xt )i | of the average total Hb copy number profile at the boundary and the boundary width ∆x with (green) and without (red) mutual repression as a function of the gap protein diffusion constant D for a system in which the gap proteins are produced in bursts of 10 at a time with decreased production rate β = β0 /10. The grey dashed lines are the values obtained from the approximation 0 ∆x = σH (xt )/|hH(xt )i |. Thick dashed lines are values for D = 0. Error bars were obtained from block averages over 10 independent samples. The black dotted line marks the D-value where the boundary is both steep and precise due to mutual repression.
23
Figure S9. The effect of increased copy number on the Hb boundary precision. See following page for description.
24
Figure S9. The effect of increased copy number on the Hb boundary precision. Shown are the value ratios of important boundary properties for production rates β > β0 as compared to β/2 for (A) the total Hb copy number noise σH (xt ) at the boundary, (B) the steepness of the average Hb profile at xt , and (C) the resulting width ∆x with (green) and without (red) mutual repression. Solid lines are for β = 2β0 , dashed lines for β = 4β0 . Blue lines depict the ratios as predicted from the expected scaling behavior for the limits of D → 0 (upper line pairs) and D −→ ∞ (lower line pairs). The steepness is expected to scale precisely with the increased copy number in both limits. Note that the expected factor of copy number increase upon doubling β is not precisely two because of the nontrivial dependence of the monomer-dimer equilibrium on the production rate.
SUPPORTING INFORMATION Mutual Repression enhances the Steepness and Precision of Gene Expression Boundaries Thomas R. Sokolowski1 , Thorsten Erdmann2 , Pieter Rein ten Wolde3,∗ 1 FOM Institute AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands 2 University of Heidelberg, Institute for Theoretical Physics, Philosophenweg 19, 69120 Heidelberg, Germany 3 FOM Institute AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands ∗ E-mail: [email protected]
1 1.1
Methodic details Measurement of the boundary width
By default we determine the boundary width in the following two ways: Let cm n,s be the copy number of species s in a nucleus with angular index m < Nφ and axial index n < Nx , where Nφ is the number of rows around the circumference of the cylinder, and Nx is the number of colums in the axial direction along the AP axis. To compute the boundary width of the expression m m = (cm domain of a gap protein s, we compute for each row m Tn,s n,s − θs ) · (cn+1,s − θs ) as a function of n, 1 m where θs is half the copy number expected at full activation. A boundary position xm t = x (nt + 2 ) is m m defined as the position (nucleus) where Tnt ,s < 0. The values of xt are recorded in a histogram; here, the positions for the different rows m are put in the same histogram. The histogram is normalized at the end of the simulation, and the boundary width ∆x is calculated as the standard deviation of this histogram. 0 Secondly, at the end of the simulation, the slope of the average, hH(xt )i , and the standard deviation of the total Hb copy number σH (xt ) at the Hb boundary position xt are calculated from the timeσ(xt ) and φ-averaged profiles. From this, an approximation for the boundary width given by ∆x ≈ |hH(x 0 t )i | is obtained, following [1–3]. To this end, first xt is determined in the same way as in the runtime measurements, only now working on the (both time- and circumference-) averaged profile. We describe 0 in the following section how the steepness hH(xt )i is measured.
1.2
Measurement of the profile steepness
In our discrete system the measurement of a local derivative at the boundary position xt is a process prone to even small stochastic variations if a naive measurement technique is chosen. If the average boundary position xt for a set of different samples with identical initial conditions always is in between two particular nuclear positions x(n0 ) and x(n0 + 1), then using linear differences to determine the 0 steepness hH(xt )i at the boundary position may give a reasonable estimate. If, however, xt fluctuates around a particular nuclear position x(n0 ) among different samples and hH(x(n0 − 1))i − hH(x(n0 ))i significantly differs from hH(x(n0 ))i − hH(x(n0 + 1))i, the linear differences method will produce a large error bar and also markedly affect the mean of xt among these samples. As a result both the measured steepness and the quality of that measurement for a given set of parameters depends on whether xt accidently happens to predominantly vary in the interval between the same nuclear positions or not. To overcome this illness we measure the boundary steepness from the average protein profile by a two-step polynomial fitting procedure: First we fit a polynomial of 3rd degree to a region of the data around xt that contains at least four points (nuclei). The derivative of the polynomial at xt gives an initial estimate of the boundary slope, which we use this to calculate the approximative x-interval over which the profile falls from maximal to minimial expression level. If the latter is larger than the original fitting range (which usually is the case) we repeat the fitting on the enlarged interval. Since the profiles to a good
2
approximation are sigmoidal functions this improves the quality of the fit. The measured boundary slope then is defined as the derivative of the polynomial function at xt after the second fitting.
1.3
Number of cortical nuclei at cell cycle 14
The development of the Drosophila embryonic syncytium starts with a single nucleus. The first 9 nuclear divisions happen in the yolk. During cell cycles 7 to 10 a migration of the nuclei towards the cortex can be observed. However, approximately 200 polyploid nuclei stay behind in the yolk and stop dividing after their 10th cycle [4]. This quiescence persists during subsequent cell cycles, including cycle 14. As an effect of this, the number of nuclei at the cortex in cycle 14 is considerably lower than 213 = 8192. An estimate of the reduced number of cortical nuclei is given by: Ncortex ' (29 − 200) · 24 = 4992
(1)
This number indeed is closer to 212 than to 213 . Note that in our model the precise number of nuclei does not matter, rather it is the distance between the nuclear compartments and the diffusion correlation length that impact on the results. Our values for both the internuclear distance and the nuclear diameter correspond to the experimental values reported by Gregor et al. [2, 5].
1.4
Predicted copy numbers and effective protein lifetime
Our main observables are the total copy numbers of Hb and Kni, defined as follows: m m H ≡ cm n,H = cn,HM + 2cn,HD + 2
5 X
cm n,K1 j
j=0 m m K ≡ cm n,K = cn,KM + 2cn,KD + 2
5 X
cm n,H1
(2)
j
j=0
Here, for G ∈ {H, K}, cm = 1 if the promoter of species G is binding j morphogen molecules and one n,G1 j
can be equal to one for only one (repressing) gap dimer; evidently, at any given moment in time cm n,G1 j
j ∈ {0..5}. The ratio between the number of monomeric and the number of dimeric proteins is a nontrivial function of the monomer production rate, the monomer and dimer degradation rates and the parameters that determine the dimerization and dedimerization reactions. To obtain an estimate for the expected copy numbers of monomers and dimers of gene g we solved the mean-field rate equations for a simplified model which comprises monomer production, (de)dimerization and monomer and dimer degradation only, i.e. in which promoter state fluctuations and diffusion are neglected, in the steady state. We assume here that stochastic monomer production events can be accounted for by an effective mean-field production rate hβi = βhH50 i for Hb and similarly for Kni, which depends on promoter (un)binding parameters and the particular morphogen and repressor levels. This yields the following prediction for the copy number
3
of monomers (GM,hβi ) and dimers (GD,hβi ): ( 1 D D 2kon µD − koff µM − µM µD GM,hβi = D 4kon µD ) q 2 D D Dµ D + 8hβikon D koff + µD + µM koff + µD − 2kon µD GD,hβi =
1 D µ2 8kon D
( D 2 D D koff µM + µD 4hβikon + µM µM − 2kon − µM
q
) D Dµ 8hβikon D koff + µD
2 D +µ D + µM koff D − 2kon µD
(3)
D D (koff ) are the dimerization forward Here µM (µD ) is the monomeric (dimeric) degradation rate and kon (backward) rates, respectively. From this we calculate the total expected copy number Ghβi := 2GD,hβi + GM,hβi at effective production rate hβi. In particular in the full-activation case, i.e. when the probability to be fully activated and unrepressed hG05 i ≈ 1 and therefore hβi ≈ β, the above estimates correspond to average values from our simulations very well. 1 (µM GM + 2µD GD ). Our stanWe define the effective degradation as µeff = µeff (GM , GD ) = GM +2G D −3 dard values result in µeff ≈ 4.34 · 10 /s with GM = GM,β and GD = GD,β .
2 2.1
Additional analysis Poissonian limit with dimerization
In [3], it was shown that, when D → ∞, the variance in the protein concentration becomes equal to the mean concentration: diffusion washes out bursts in gene expression, thus reducing the non-Poissonian part of the noise. However, in that model the proteins do not dimerize, in contrast to our model. With dimerization, a different limit for the variance in the total protein concentration is approached as D → ∞. To derive this limit, first note that the total protein copy number G of a protein G is G ≈ 2GD + GM . Assuming that hGD GM i ≈ hGD ihGM i (our simulations indicate that this approximation is very accurate), 2 we find that the variance σG in G is: 2 2 2 σG ≈ 4σG + σG , D M
(4)
2 2 where σG is the variance in the dimer level GD and σG is the variance in the monomer level GM . Both D M monomers and dimers are subject to spatial averaging, and therefore their variances can be written in the form [3]:
1 2 σ0,G − GM M N 1 2 = GD + σ0,G − GD D N
2 σG = GM + M 2 σG D
(5)
2 Here N is the number of nuclei contributing to the averaging, which is proportional to D, and σ0,G is M/D the variance in the monomer and dimer levels in the absence of diffusion, respectively. The part preceded by 1/N represents the variance that can be reduced by spatial averaging. Plugging these expressions into
4
the previous and using G = 2GD + GM we arrive at: 1 2 2 2 4σ0,GD + σ0,G − 4GD − GM σG = 4GD + GM + M N 2GD 1 2GD 2 = 1+ G+ σ0,G − 1+ G G N G 1 2 σ =: (1 + fD ) G + − (1 + fD ) G N 0,G
(6)
Note that N is the same for both monomers and dimers because their diffusion√constant does not differ in our p model. Evidently, the lower bound for σG in the limit N → ∞ is not G any more, but given by (1 + fD ) G, where fD is the fraction of proteins in the dimer state with respect to the total protein number (implying fD ≤ 1). This is indeed what we observe in our data for σG . In our simulations the equilibrium is strongly shifted towards the dimerized state, so that fD ≈ 0.97. We can understand the limit N, D → ∞ intuitively by noting that in this limit there is no noise in the nuclear protein concentration due to the stochastic production and decay of molecules in each of the nuclei—this is because the synthesized molecules are immediately donated to a reservoir that is infinitely large; instead, there is only noise in the nuclear protein concentration due to the sampling of molecules from this 2 2 reservoir, which obeys Poissonian statistics: σG = GM and σG = GD . This yields, for N, D → ∞, M D 2 2 2 σG = 4σGD + σGM = 4GD + GM = (1 + fD )G.
2.2
Bifurcation analysis
In order to predict the regions in which bistability can be expected for different amplitudes A of the morphogen gradients we performed a deterministic mean-field bifurcation analysis for a simplified 1dimensional version of our model of mutual repression between hb and kni . The analysis is based on the following two equations describing the change of the mean-field total copy number of Hb (H(x)) and Kni (K(x)) at position x: ∂t H(x) = βH (x) ∂t K(x) = βK (x)
2 KR 2
− µH H(x)
(7)
2
− µK K(x)
(8)
2 + [f K(x)] KR D 2 KR
2 + [f H(x)] KR D
Here βH and βK represent the protein synthesis rates, µH and µK the corresponding (effective) degradation rates, KR is the dissociation constant of cooperative repressor binding to the promoter and fD is the fraction of proteins in the dimerized qstate. Note that, since the intermediate step of dimerization is
R /k R if k R and k R are the binding rates of the dimers. To neglected here, we have to take KR = koff on on off facilitate calculations we make two further simplifying assumptions here:
1. We neglect activation dynamics and resulting promoter state fluctuations, i.e. we assume that certain constant levels of the activators at position x lead to average constant production rates βH = β([Bcd](x)) and βK = β([Cad](x)), respectively. In our standard case β([Act](x)) = [Act]5 (x)/([Act]5 (x) + 6905 ) for both [Act] = [Bcd] and [Act] = [Cad]. 2. In our simulations we have different degradation rates for monomers and dimers so that the effective total degradation rate depends on the monomer-to-dimer ratio, which in turn varies with the total copy number (see section 1.4). Thus, in principle, also fD and µH and µK are functions of x, or the corresponding activator levels. Since this introduces further nonlinearities into the above equations and complicates their solution, we substitute the degradation rates µH and µK by a
5
constant value µeff , which is the effective degradation rate for the maximal expression level (full activation). Also for fD we take the constant value for full activation, fD ' 0.97, which reflects that the dimerization equilibrium in our simulations is strongly shifted towards the dimerized state. The predictions concerning the bifurcation behavior only change marginally if µeff and fD values for lower expression levels are used. For each position x with local activator levels corresponding to the ones in the simulations we calculated fixed point solutions for the copy number pair (H(x), K(x)) starting from the steady-state assumption ∂t (H(x), K(x)) = (0, 0). The stability of the fixed points was determined starting from the Jacobian for the above ODE system: ∂H [∂t H] ∂K [∂t H] J(H, K) = (9) ∂H [∂t K] ∂K [∂t K] Within the relevant parameter regime we obtained fixed points with either two negative eigenvalues (i.e. stable fixed points) or one positive and one negative eigenvalue (i.e. saddle points). The determinant therefore completely characterizes the stability of the fixed points. If det J(H0 , K0 ) < 0, then (H0 , K0 ) is a saddle point. Otherwise it is stable. Fig. S1 shows the fixed point solutions for Hb and Kni as a function of x for different activator amplitudes A. Stable solutions are drawn with solid, unstable solutions with dashed lines. Depending on the A value, the system displays a saddle node bifurcation at a point towards the anterior (Hb) or posterior (Kni) from midembryo. Within the region confined by the bifurcation points two stable and one unstable fixed points exist for each gene, implying bistability. The region clearly widens for increasing A and spans almost the whole embryo length for A = 8. Our deterministic analysis therefore predicts the enlargement of the region of bistability as observed in our single nucleus simulations.
2.3
Estimation of switching times
To quantify the swiching times in the presence of bistability we performed simulations of isolated single nuclei featuring the same set of reactions and parameters as in the full scale simulation. To obtain estimates of switching times at different positions x along the AP axis we set the levels of Bcd and Cad in the given nucleus equal to the ones at x in the space-resolved simulations. The switching time was estimated by calculating from long time trajectories of the total Hb and Kni copy numbers the relaxation time ts of the average correlation function hC(t)it0 ≡
hIH (t0 )IK (t)it0 hIH (t0 )it0
(10)
where IH (IK ) are indicator functions which are one if the difference in the total gap gene copy numbers ∆N = H − K is above (below) a certain threshold ΘN (−ΘN ). ΘN thus defines the regions within which the switch is considered to have switched to the Hb–high or Kni–high states, respectively, and serves to separate the stable attractor states from the transition region. We found that ΘN = 200 is a reasonable choice for our set of parameters. We determined the switching times from one long sample for different positions x and different activator amplitudes A and find that ts is very similar within the double-activated bistable regions for high A. To obtain an error estimate we additionally calculated block averages of estimated switching times among 10 long samples for various A at midembryo (x = L/2). Table S1 shows our results from the latter procedure. Note that for A = 1 the system is not truly bistable yet because for A = 1 we have half-activation at midembryo and due to the lack of diffusion large promoter-state fluctuations dominate over long-time switching potentially induced by mutual repression. Consequently, the given number does not reflect a switching time. We cite it here for completeness, however.
6
Table S1. Switching times at midembryo for different activator levels. Activator amplitude A Switching time ts [s] (1) (6343.7 ± 17.2) 2 20302.5 ± 74.5 4 20957.3 ± 54.9 8 20994.7 ± 67.2
2.4
Analysis of statistical properties of the boundary
Our measures for both the boundary steepness and the variance of the boundary position are based on averages over both the time and the circumference of the embryo which were calculated during runtime. While the double-averaging procedure limits the amount of data that must be stored and facilitates rapid acquisition of good statistics, it also discards information about the microscopic properties of the boundary at a given time instance. Based on the average data it is impossible to determine whether the blurring of the boundary quantified by ∆x is due to concerted stochastic movements of a steep and rather homogeneous instantaneous boundary or simply due to stochastic fluctuations of the boundary position in each nuclear row around a well-defined constant mean boundary position (or due to both). In the latter case the boundary will be rough at each given time instance, i.e. the time average of the boundary position variance in the cicumferential direction will be large, but the time variance of its circumferential mean will be negligible. The opposite will be the case in the other extreme. These quantities therefore can be used to distinguish the two hypothetical situations. The overall boundary width in both cases is given by the sum: 2 ∆x2 = σx2t (φ) + σhx t iφ (t)
(11)
Here h. . . iφ denotes the average over the circumference, while the bar denotes the time average. An identical variance decomposition can be made for the fluctuations of the Hb copy number at any position x along the AP axis. Similarly, comparing the average of the profile steepness for a particular nuclear row and time instance to the steepness of the time- and circumference average of the copy number reveals whether the steepness of the average profile is due to concerted movements of similarly shallow instantaneous profiles or due to unconcerted fluctuations of steep instantaneous profiles. In order to determine which of the portrayed blurring mechanisms is dominant in our system we performed the described variance decomposition for a set of 100 instantaneous outputs of the fully resolved 2D system in steady state, i.e. for 6400 different total Hb copy number profiles along the AP axis, for both the variances at the boundary and for the steepness at the boundary and for both the system with and without mutual repression. We focused on our standard parameter set (see Table S2) and a range of gap protein diffusion constants D. 2.4.1
At each time instance the boundary is rather rough
Fig. S2 shows for the systems with (Fig. S2A and C) and without (Fig. S2B and D) mutual repression the variance decomposition for the variance of the Hb copy number at the boundary (Fig. S2A and B) and for the variance of the boundary position xt (Fig. S2C and D) as a function of the Hb protein diffusion constant D. As a control we compare the total variances calculated from the instantaneous profiles to the variances accumulated during runtime and, in case of ∆x, to the value obtained from the 0 approximation ∆x = σH (xt )/|hH(xt )i | (note that here h. . . i is the average over both time and φ). We see a good agreement between these quantities. The plots reveal that both for σH (xt ) and ∆x the variance over the circumference at a fixed time is by far the dominant contribution to the overall variance. This implies that in our system the boundary is indeed very rough at each time point and that concerted boundary movements do not occur.
7
2.4.2
At each time instance the profiles are slightly steeper than their average
The calculation of the variance decomposition is less straightforward for the slope. In particular for low D, when spatial averaging is still inefficient, the instantaneous profiles are very ragged and the boundary threshold value typically is crossed at multiple positions along the AP axis. This makes it impossible to uniquely define an instantaneous boundary position as required to calculate the instantaneous boundary slope. In order to perform the analysis at least on a subset of the data we introduced a protocol which only takes into account instantaneous profiles with a single boundary crossing, rejecting all other profiles. For low D, however, the rejection rates rise above 90%. We therefore decided to smoothen the profiles by computing running averages between a fixed number ν of nuclei along the AP axis before the analysis. The averaging lowers the rejection rate dramatically, however it also decreases the profile steepness and therefore manipulates the observable of interest. Nevertheless we can make a qualitative statement on the base of the results obtained for only slight smoothening of the profiles (ν = 3). For simplicity and due to increased data abundance, in this analysis we used simple finite differences to determine the slope. In Fig. S3 we plot the average of the instantaneous boundary steepness for different degrees of smoothening (averaging over ν = 3, 5, 7 nuclei along the x-axis) as a function of D and compare to the steepness of the average Hb profile for the system with (S3A) and without (S3B) mutual repression. While the data for ν > 3 clearly must be considered biased by the running averages, the values for ν = 3 show that the instantaneous boundary slope on average is higher than the slope of the average profile, in particular for low diffusion constants. The variance decomposition for the boundary position xt shows that the variance of the circumference mean of the boundary position in time is very small. This implies that the steepness of the circumferenceaveraged profiles should be approximately equal to the steepness of the time- and circumference-averaged profile. As a control we therefore repeated the above analysis on the 100 φ-averaged instantaneous profiles of the same dataset. The averaging along the circumference significantly reduces the number of profiles with ambiguous boundary positions. We therefore were able to obtain reasonable estimates of the observable without pre-smoothening of the profiles (ν = 1). The results are shown in Fig. S3C for the system with mutual repression and Fig. S3D for the system without mutual repression. In the system with mutual repression the average slope of the φ-averaged profiles for ν = 1 agrees well with the slope of the both time- and φ-averaged Hb profile. In the system without mutual repression the φ-averaged profiles are slightly steeper than the average.
3
Additional simulations
3.1
Influence of the Hill coefficient
To address the influence of changing activator cooperativity on our results we performed simulations with reduced number of activator binding sites nmax . While in our model this is achieved by simply reducing the number of intermediate states between the empty promoter state and the producing promoter state, the binding parameters have to be rescaled with care to preserve the activation equilibrium at midembryo. Since we assume the activator binding rates to be diffusion limited, the necessary changes affect the A unbinding rates koff,n = a/bn . However, even when preserving the equilibrium, the freedom in the choice of these parameters allows for altering the time scale of transitions between the different activation levels. In order to rescale the rates in a unique fashion upon lowering nmax we imposed the following constraints: A 1. For all nmax the effective activator dissociation constant at midembryo KD = A1/2 = 690 is preserved, which implies that for all nmax the average activation probability at midembryo is 1/2.
2. The waiting-time distribution for the unbinding from the producing state is the same for all nmax A and, for comparison, equal to the one for the default cooperativity nmax = 5, i.e. ∀n : koff,n = max A const = koff,5 .
8
3. The off-rate reduction per subsequent activator binding is always 1/b, i.e. ∀n : b(n) = b. A A A Note that for nmax = 1 the first two conditions can be met together only if KD kon = koff,5 , which is not the case for our parameter set. We therefore restricted ourselves to nmax ∈ {2, 3, 4, 5}. For each nmax , the above constraints were used to uniquely determine the parameters a and b from the exact analytical solution for the average occupancy of the producing state, which was obtained from steadystate mean-field solutions of the chemical mass-action ODEs. Interestingly this results in only minor differences in a among the different nmax values, while the reduction per binding step 1/bnmax becomes significantly larger for lower nmax . This fact has an important implication for the noise charactistics of the different promoters: If anmax ' a = const for all nmax then the unbinding rate from the state binding (nmax − 1) activator proteins (the “highest” non-producing activator state) is given by: A A max −1) koff,(n ' a/bn(nmax = bnmax koff,5 max −1) A Since bnmax markedly increases with decreasing Hill coefficient the off-rate koff,(n for low nmax will max −1) be higher than the corresponding rate for high nmax . This will favor rapid returns to the producing state with nmax bound activator molecules for high Hill coefficients, whereas for low Hill coefficients the promoter is more likely to descent into the regime with less activator molecules bound. The fact that this is less likely for higher Hill coefficients is compensated by the fact that also the time to return to the producing state from the states binding low numbers of activator molecules on average is longer for higher nmax . Note that the mean off-time–just as the mean on-time–is the same for all nmax . In short, for the promoters with higher Hill coefficients we expect an off-time distribution with high probability weight on short off-times and a long low-probability tail for long off-times, while the distribution for lower Hill coefficients should resemble an exponential. In order to illustrate this effect we recorded long time-trajectories of the occupancy of the producing state in a single isolated nucleus close to midembryo for different nmax and without mutual repression nor diffusion. All other parameters were kept at the standard values. From these trajectories we determined the on- and off-times of the promoter and binned them into a histogram. The results are shown in Fig. S5. It can be seen that while for nmax = 2 the two distributions are exponential with approximately equal mean, the off-times distribution increasingly deviates from an exponential distribution as nmax is increased; more probability is shifted to very short off-times and very long off-times, causing the emergence of a long tail in the distribution.
3.1.1
Also for lower Hill coefficients mutual repression steepens profiles without corrupting boundary precision
The broadening of the off-times distribution is expected to result in higher output noise for high nmax as compared to low nmax . This is confirmed by the simulations of the full-scale spatially resolved system for different nmax . Fig. S6 shows σH (xt ), the average standard deviation of the total Hb copy number at 0 the boundary position xt (upper panels), the steepness |hH(xt )i | of the average Hb profile at xt (middle panels) and the boundary width ∆x (lower panels) as a function of the gap protein diffusion constant D for nmax ∈ {2, 3, 4, 5}. σH (xt ) is indeed decreasing upon lowering nmax , in particular in the regime of low diffusion constants. For higher diffusion constants the decrease is less pronounced: spatial averaging is efficient enough to lower the output noise down to the observed values irrespective of the width of the offtime distribution. The noise decreases less markedly for the systems with mutual repression. This is most likely due to the fact that lowering nmax also increases the probability of occasional repressor production beyond midembryo, which in turn increases the noise. The steepness plots reveal that, although the profiles naturally become less steep for lower nmax , the steepness in the systems with mutual repression is markedly higher than the one in the system without mutual repression. In all systems the steepness as a function of D shows a very similar behavior: Upon increasing D the steepness in the systems with mutual repression first increases towards a maximum before it rapidly decreases. Since both σH (xt ) and
9 0
|hH(xt )i | change with nmax in a similar fashion, in particular in the region around D = 1 µm2 /s, the width ∆x as a function of D also looks very similar in this region for all nmax . In all cases the profiles in the system with mutual repression are more precise and markedly steeper as compared to the system without mutual repression at a D-value which is one order of magnitude less than the optimal value in the systems without mutual repression. Therefore the basic effect observed in our simulations for nmax = 5 persists in the simulations for lower Hill coefficients. 3.1.2
Lower Hill coefficients allow for stronger morphogen level variations
Although lowering nmax in our system reduces the protein production noise it also markedly decreases the steepness of the gene activation profiles. An important implication of this is that for lower nmax the activation probability beyond midembryo increases. Lowering nmax thus is similar to increasing the activator amplitude A and, in principle, might result in the creation of a bistable region around midembryo already for lower A-values as compared to the system with nmax = 5. We analysed how the results for ∆x as a function of A for the case of correlated variations change as nmax is decreased. Fig. S7 shows ∆x(A) for nmax ∈ {2, 3, 4, 5} and D = 1.0 µm2 /s for systems with and without mutual repression. Overall, ∆x(A) is very similar for all considered nmax . For A ≤ 2 the width ∆x in the systems with mutual repression is always lower than in the systems without mutual repression. The minimal ∆x is attained at A = 1 in all cases. The main difference is in how ∆x changes with A for A > 1: The lower nmax , the slower the width increases with A. Thus, while lower Hill coefficients decrease the steepness of the profiles significantly, they may prove beneficial by extending the range over which extrinsic variations are successfully buffered without increasing intrinsic fluctuations of the boundary.
3.2
Influence of the expression level
In order to examine the influence of a changed signal-to-noise ratio on our results we performed simulations with altered production dynamics. We did this by (1) introducing bursty production, i.e. producing 10 copies of the gap protein monomer at a time with a 10 times lower production rate (β = β0 /10), and (2) by keeping the burst size at one and changing the production rate β. To preserve the binding equilibrium of the repression reaction at midembryo upon changing β we also changed the off-rate of the repressor dimers by a factor fβD , which is the ratio between the expected number of dimers at midembryo for the altered production rate β and the corresponding value for the standard production rate β0 . Note that, since in our system the copy numbers of both monomers and dimers depend on β in a nontrivial fashion (see section 1.4) the effective copy number increase typically does not correspond to the ratio β/β0 . Therefore fβD > β/β0 for β > β0 . 3.2.1
Bursty production has only a marginal influence on the boundary properties
In Fig. S8 we plot the standard deviation of the total Hb copy number at the boundary, the steepness of the total Hb copy number profile at the boundary and the boundary width ∆x as a function of D for the system with bursty production (burst size 10). There is no significant difference as compared to the system with normal production (burst size 1, compare to Fig. S6(D) or Fig. 3 of main article). For low D the production noise is marginally higher with bursty production, resulting in a slight increase of ∆x in this regime; the effect of varying D, however, is much more important. This is most likely a consequence of the fact that for the given Hill coefficient nmax = 5 promoter state fluctuations are already at a high level due to a very broad off-time distribution (see section 3.1).
10
3.2.2
Increased production rates reveal different noise scaling behavior for different regimes of the diffusion constant
Upon increasing the production rate and consequently the total copy number of the gap proteins we may expect a relative decrease in the output noise, but only if the latter is purely Poissonian. In our system √ this corresponds to the limit of high gap protein diffusion constants. In that limit, we expect σG ∝ G, where σG is the noise in the total gap gene copy number G. However, in the abscence of spatial averaging, i.e. for the limit D → 0, non-Poissonian noise prevails and the expected scaling is σG ∝ G [3]. If the copy number profile is scaled uniformly at each AP position x, which–to a good approximation–is the case in our system, we expect for the scaling of the gradient at midembryo G0 (xt ) ∝ G(xt ). The expected √ scaling for the boundary width ∆x then is ∆x ∝ 1 for low diffusion constants and ∆x ∝ 1/ G for high D. While the overall characteristics of the boundary are very similar to the system with β = β0 , a comparison roughly confirms the predicted scaling. Fig. S9 compares for Hb the standard deviation of the total copy number at the boundary (S9A), the steepness at the boundary (S9B) and the resulting boundary width (S9C) as a function of D for increased production rates β = 2β0 and β = 4β0 to the corresponding values for the sytem with production rate β/2. Thus, the values for β = 4β0 are compared to β = 2β0 and the values for β = 2β0 are compared to β0 . Blue lines mark the expected change of the quantities as predicted by the scaling relations, where the corresponding copy number increase is given by the factor f2 ≡ f2β0 = 2.22 and f4 ≡ f4β0 /f2β0 = 2.16. Here fβ ' fβD is the predicted total copy number at midembryo divided by the corresponding value for β = β0 . The plots show that while the slope ratio is roughly equal to f2 (f4 ) for all D both in the system with (green) and without (red) mutual repression, the noise ratio depends on the diffusion constant and also slightly differs for the systems with and without mutual repression. Nevertheless the predicted scaling behavior is confirmed√in both √ cases: in the low diffusion constant regime the noise ratio is roughly f2 together this leads to a boundary width ratio of one for (f4 ) and approaches f2 ( f4 ) as√D increases; √ low D which decreases towards 1/ f2 (1/ f4 ) for higher D.
3.3
Activation of both gap genes by a single gradient
In the one-morphogen gradient scenario, both hb and kni are activated by the Bcd gradient. Here, kni is activated in the same way as hb, namely by 5-step cooperative binding, but with a lower activation threshold. This results in the induction of both genes in the anterior half of the embryo up to the posterior Hb boundary and of kni in an additional region posterior to the Hb boundary. Given that hb represses kni more strongly than vice versa in the double-activated bistable region, this parameter choice will result in the formation of two neighboring domains. We chose the kni activation threshold to be lower by a factor of 1/2, which causes an offset of its half-activation point by approximately 10 nuclei (83 µm) towards R,K the posterior. We varied the protein diffusion coefficient D and koff , the off-rate of the Kni repressor dimers from the hb promoter. The rate for dissociation of the Hb dimers from the kni promoter was kept R,H at the standard value koff = 5.27 · 10−3 /s in all simulations. R,K The diffusion constant D of the gap proteins and the dissociation rate koff of Kni from the hb promoter are indeed key parameters. On the one hand, hb must repress kni more strongly than the other R,K way around, because otherwise there will be only one kni domain. On the other hand, when koff is high, then kni is only significantly expressed when D is low, because kni represses hb more weakly than vice versa, which means that low amounts of invading Hb dimers are sufficient to shut off Kni production almost completely; indeed, in this regime, kni has hardly any effect on the precision of the hb expression R,H R,K domain. We found that when koff /koff is roughly between 0.1 and 1, both hb and kni domains are √ R,H R,K formed robustly. In Fig. S4 we display the case for koff /koff = 1/ 10 ≈ 1/3. Fig. S4A shows that the maximum of the average Kni copy number is lower than that of Hb, even though for x < 60 %EL kni is essentially fully activated by Bcd (Fig. S4B). The lower maximum is due to the fact that hb represses kni more strongly than vice versa. Another point worthy of note is that the
11
fluctuations in the Kni copy number in the Kni domain are higher than those of Hb in the Hb domain (Fig. S4C). This is essentially due to the small width of the Kni domain: kni is either fully activated by Bcd yet still repressed by Hb or not repressed by Hb yet stochastically activated by Bcd. Panels D-F show, respectively, the noise in the Hb copy number at the hb expression boundary, the steepness of this boundary, and the width of this boundary, as a function of the diffusion constant D of the gap proteins. It is seen that the results are highly similar to those of the two-gradient motif. The noise in the Hb copy number at the boundary is not much affected by mutual repression (panel D), while the steepness, and consequently boundary precision, is markedly enhanced by mutual repression, especially when the diffusion constant is small. Note that while for the two-morphogen gradient scenario 0 the approximation ∆x ≈ σH (xt )/|hH(xt )i | is in very reasonable agreement with ∆x as measured from the distribution of threshold crossings p(x), here the agreement is much less. This is due to sporadic repression events in the anterior region where hb and kni are both fully activated, which leads to a long tail of p(x) extending towards the anterior pole; while p(x) in the tail is small, the fact that the tail is long does 0 markedly increase the standard deviation ∆x. Given that the approximation ∆x ≈ σH (xt )/|hH(xt )i | works so well for all the other cases, we consider this approximation, which does not suffer from sporadic but strong hb repression events in the anterior, to be more reliable. We therefore conclude that also in the one-morphogen gradient scenario, mutual repression can enhance both the steepness and the precision of gene-expression boundaries.
References 1. Tostevin F, ten Wolde PR, Howard M (2007) Fundamental limits to position determination by concentration gradients. PLoS Comput Biol 3: e78. 2. Gregor T, Tank DW, Wieschaus EF, Bialek W (2007) Probing the limits to positional information. Cell 130: 153–164. 3. Erdmann T, Howard M, ten Wolde PR (2009) Role of spatial averaging in the precision of gene expression patterns. Phys Rev Lett 103: 258101. 4. Foe VE, Alberts BM (1983) Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in drosophila embryogenesis. J Cell Sci 61: 31–70. 5. Gregor T, Wieschaus EF, McGregor AP, Bialek W, Tank DW (2007) Stability and nuclear dynamics of the bicoid morphogen gradient. Cell 130: 141–152. 6. Abu-Arish A, Porcher A, Czerwonka A, Dostatni N, Fradin C (2010) High mobility of bicoid captured by fluorescence correlation spectroscopy: Implication for the rapid establishment of its gradient. Biophys J 99: L33–L35.
Table S2. The standard set of the Parameter / Quantity Symbol Geometry Number of nuclei in axial direction x Nx Number of nuclei in circumferential direction φ Nφ Internuclear distance (lattice spacing) l - resulting embryo length L Nuclear radius R - resulting nuclear volume V Morphogen / activator gradients Standard copy number at the poles A0 - resulting copy number at half activation A 12 Morphogen diffusion length λ Activation dynamics Intranuclear activator diffusion constant DA Activator binding site target size α A - resulting activator binding rate kon A Activator unbinding rate koff,n Hill coefficient of activation nmax Gap gene dynamics Gap protein monomer production rate β0 Gap protein monomer degradation rate µM Gap protein dimer degradation rate µD D Dimerization forward rate kon D Dimerization backward rate koff - predicted total Hb (Kni) copy number at full activation and D = 0 Hβ0 (Kβ0 ) Gap protein diffusion constant D Repression dynamics R Repressing dimer binding rate kon R Repressing dimer unbinding rate koff 0.40 5.27·10−3
µm3 /s 1/s
1/s 1/s 1/s µm3 /s 1/s
3.37 3.37·10−2 3.37·10−3 0.80 5.59·10−3 775 varied
µm2 /s nm µm3 /s 1/s
µm
6720 690 119.5 3.2 10 0.40 410/6n 5
µm µm µm µm3
64 64 8.5 544 6.5 143.8
equal to 0.1 ·
A same as kon varied in some simulations − kA,5 ;
A same as 2kon
A = 4παDA diffusion limited, kon n = no. of bound activator molecules supported by [2]
supported by [6]
varied in some simulations supported by [2] supported by [2]
supported by [5] assuming spherical shape
supported by [2]
most important parameters in our model. Value Unit Remarks
12
13
Figure S1. Bifurcation analysis of a one dimensional mean-field model of mutually repressing gap genes activated by morphogen gradients. Plotted are the stable (solid lines) and unstable (dashed lines) fixed points of the copy number of Hb (colored lines) and Kni (grey lines) as a function of the AP position x as predicted by a bifurcation analysis performed on a 1D mean-field model in which hb and kni are activated cooperatively by their respective morphogens and mutually repressing each other. Different colors correspond to different fixed points. The different panels show the solutions for activator amplitudes (A) A=1, (B) A=2, (C) A=4 and (D) A=8. All other parameters values are the standard values from Table S2. Activator concentrations at position x used in the mean-field analysis correspond to the ones in the 2D stochastic simulations. Away from midembryo each gap protein level has only one stable fixed point and one of the two levels is always zero. For all A there is a region around midembryo in which the protein levels have two stable and one unstable fixed point, implying bistability. In this region the analysis predicts bistable switching between the high-Hb–low-Kni and the low-Hb–high-Kni state. For clarity we color-code the Hb fixed points only. The Kni solutions are identical to the Hb solutions mirrored with respect to midembryo.
14
Figure S2. Decomposition of variances at the boundary. (A) Decomposition of the total Hb copy number variance at the average boundary position for the system with mutual repression as a 2 function of the gap protein diffusion constant D. Plotted are: σH (xt ) the total (time- and circumference-) variance measured during runtime (RT, black), the same quantitity determined from a 2 set of 6400 instantaneous profiles (IP, red, 64 AP rows at 100 different time points), σhHi the variance φ
2 the time average of the variance of in time of the circumference average of H(xt , φ) (green) and σH H(xt , φ) over the circumference (blue). (B) The same as (A) for the system without mutual repression. (C) The same variance decomposition as in (A) for the Hb boundary position xt instead of the copy number. The black line shows the ∆x values measured as the standard deviation of the boundary position histogram accumulated during runtime (RT), the grey dashed line the corresponding values 0 determined from the approximation σH (xt )/|hH(xt )i |. (D) The same as in (C) for the system without mutual repression. In both cases, the main contribution to the total boundary variance σx2t comes from σx2t , implying that the blurring of the boundary is rather due to roughness than due to concerted boundary movements.
15
Figure S3. Microscopic properties of the boundary steepness. Panels (A) and (B) compare for different gap protein diffusion constants D the average Hb profile steepness at the boundary measured in a set of 6400 instantaneous profiles (64 AP rows at 100 different time points) to the 0 steepness of the (time- and circumference-) average of the Hb profile (|hH(xt )i |, black) for different numbers ν of neighboring data points used in calculating running averages over the instantaneous profiles for the system with (A) and without (B) mutual repression. Although for increasing ν the instantaneous profiles become less steep as a consequence of the smoothening, the values for ν = 3 indicate that the profiles at a given row and time instance are slightly steeper than the average profile. In panels (C) and (D) we show results of the same analysis performed on the 100 circumference-averages of the instantaneous profiles, again for the system with (C) and without (D) mutual repression. Here ν = 1 is the data obtained without calculating running averages (magenta). In both systems the steepness of the φ-averaged profiles agrees reasonably well with the steepness of the 0 average profile |hH(xt )i |.
16
Figure S4. The effect of mutual repression in a system where both hb and kni are activated by the Bcd gradient. See following page for description.
17
Figure S4. The effect of mutual repression in a system where both hb and kni are activated by the Bcd gradient. (A) Time- and circumference-averaged Hb (hHi, solid lines) and Kni (hKi, dashed lines) total copy-number profiles along the AP-axis for various Hb and Kni diffusion constants D. (B) AP profiles of the average standard deviation of the total gap gene copy number for Hb (σH , solid lines) and Kni (σK , dashed lines). Note that the noise in K in the Kni domain is larger than that in H in the Hb domain. (C) AP profiles of the probabilities hH50 i + hH51 i and hK50 i + hK51 i that the hb (solid lines) and kni (dashed lines) promoters have 5 Bcd molecules bound to them, respectively; in the absence of mutual repression between hb and kni , these profiles would directly determine the expression of hb and kni . (D) The noise in the Hb copy number at the hb expression boundary as a function of the Hb and Kni diffusion constant D. (E) The steepness of the hb expression boundary as a function of the diffusion constant of the gap proteins. (F) The width ∆x of the hb expression boundary as a function of the Hb and Kni diffusion constant. The grey line corresponds to 0 the approximation ∆x ≈ σH (xt )/|hH(xt )i |, which we consider to be more reliable than ∆x as measured from the distribution of threshold crossings, p(x); the latter suffers from sporadic but strong suppression events of hb by kni in the anterior, which leads to a long tail of p(x), increasing ∆x. It is seen that while mutual repression has hardly any effect on the noise in the copy number at the boundary, it does markedly enhance the steepness of the boundary, and thereby its precision. The ratio of the R,H R,K Hb–kni -promoter dissociation rate over the Kni–hb-promoter dissociation rate is koff /koff = 1/3.
18
Figure S5. On- and off-times distributions of the hb promoter for different Hill coefficients nmax in a nucleus at midembryo. The panels show normalized histograms of the times spent in the producing (n = nmax ) promoter state (“ON”, green) and of the times spent in the non-producing (n < nmax ) states (“OFF”, red) for (A) nmax = 2, (B) nmax = 3, (C) nmax = 4 and (D) nmax = 5 (standard case). It can be seen that with increasing Hill coefficient nmax the off-times distribution changes from an exponential to a non-exponential distribution with high weight on very short off-times (implying fast returns to the producing state) and a with a long tail of long off-times. Since the off-rate from the producing state is kept the same for all nmax the on-times distributions remain unaltered. The on- and off-times have been determined from long time trajectories (ttotal = 105 s) of the occupancy of the producing state with a sampling resolution of 0.5 s.
19
Figure S6. Boundary characteristics for reduced Hill coefficients nmax . See following page for description.
20
Figure S6. Boundary characteristics for reduced Hill coefficients nmax . The standard deviation of the total Hb copy number at the boundary (σH (xt ), upper panels), the gradient of the 0 average Hb total copy number gradient at the boundary (|hH(xt )i |, middle panels) and the boundary width (∆x, lower panels) as a function of the gap protein diffusion constant D for the systems with (green) and without (red) mutual repression and Hill coefficients (A) nmax = 2, (B) nmax = 3, (C) nmax = 4 and (D) nmax = 5 (standard case). Grey dashed lines are values determined from the 0 approximation ∆x = σH (xt )/|hH(xt )i |, solid lines are values calculated from the distributions of xt . Broad dashed lines are the values for D = 0. Black dotted lines mark the D-value where the boundaries are both steep and precise due to mutual repression.
21
Figure S7. The effect of changing the activator amplitude A on the boundary precision for reduced Hill coefficients nmax . Shown are the the boundary width ∆x with (green) and without (red) mutual repression as a function of ∆xA , the separation between the Hb and Kni boundaries expected in the system without mutual repression, and the corresponding activator amplitude A for Hill coefficients (A) nmax = 2, (B) nmax = 3, (C) nmax = 4 and (D) nmax = 5 (standard case). In all cases D = 1.0 µm2 /s.
22
Figure S8. The effect of bursty gap protein production on the Hb boundary precision. The plot shows σH (xt ) the standard deviation of the total Hb copy number at the boundary, the 0 steepness |hH(xt )i | of the average total Hb copy number profile at the boundary and the boundary width ∆x with (green) and without (red) mutual repression as a function of the gap protein diffusion constant D for a system in which the gap proteins are produced in bursts of 10 at a time with decreased production rate β = β0 /10. The grey dashed lines are the values obtained from the approximation 0 ∆x = σH (xt )/|hH(xt )i |. Thick dashed lines are values for D = 0. Error bars were obtained from block averages over 10 independent samples. The black dotted line marks the D-value where the boundary is both steep and precise due to mutual repression.
23
Figure S9. The effect of increased copy number on the Hb boundary precision. See following page for description.
24
Figure S9. The effect of increased copy number on the Hb boundary precision. Shown are the value ratios of important boundary properties for production rates β > β0 as compared to β/2 for (A) the total Hb copy number noise σH (xt ) at the boundary, (B) the steepness of the average Hb profile at xt , and (C) the resulting width ∆x with (green) and without (red) mutual repression. Solid lines are for β = 2β0 , dashed lines for β = 4β0 . Blue lines depict the ratios as predicted from the expected scaling behavior for the limits of D → 0 (upper line pairs) and D −→ ∞ (lower line pairs). The steepness is expected to scale precisely with the increased copy number in both limits. Note that the expected factor of copy number increase upon doubling β is not precisely two because of the nontrivial dependence of the monomer-dimer equilibrium on the production rate.
