 Article
 Figures and data
 Abstract
 eLife digest
 Introduction
 Analysis
 Discussion
 Materials and methods
 Appendix 1: Detailed balance
 Appendix 2: Recurrence example
 Appendix 3: Application of recurrence relation to {de novo; decay} model
 Appendix 4: Application of recurrence relation to {fission; fusion} model
 Appendix 5: Application of recurrence relation to {de novo, fission; decay} model
 References
 Decision letter
 Author response
 Article and author information
 Metrics
Abstract
We present a reanalysis of the stochastic model of organelle production and show that the equilibrium distributions for the organelle numbers predicted by this model can be readily calculated in three different scenarios. These three distributions can be identified as standard distributions, and the corresponding exact formulae for their mean and variance can therefore be used in further analysis. This removes the need to rely on stochastic simulations or approximate formulae (derived using the fluctuation dissipation theorem). These calculations allow for further analysis of the predictions of the model. On the basis of this we question the extent to which the model can be used to conclude that peroxisome biogenesis is dominated by de novo production when Saccharomyces cerevisiae cells are grown on glucose medium.
https://doi.org/10.7554/eLife.10167.001eLife digest
Any cell that has a nucleus also contains a number of subcellular structures called organelles. The number of organelles inside a cell increases when new organelles are made from scratch (a process known as de novo synthesis), or when an existing organelle divides to produce two organelles in a process called fission. And the number of organelles decreases when an existing organelle decays, or when two organelles fuse together to become one organelle. The actual number of organelles of a particular type inside a cell results from a balance between these creative and destructive processes.
Last year researchers at Harvard University developed a model that treats the processes of organelle creation and destruction as if they were chemical reactions, and then used their model to make predictions about the budding yeast S. cerevisiae in three scenarios. The Harvard researchers had to use a number of approximations to make these predictions.
Now Jeremy Craven has derived exact solutions to the model for these three scenarios. The exact solutions call into question some aspects of the model, notably the prediction that the production of new peroxisomes – organelles that are involved in breaking down fatty acids and other compounds – is dominated by fission when the yeast cells are grown on a substance called oleate, and by de novo synthesis when they are grown on glucose. Craven's analysis also highlights the need for quantitative timecourse imaging data to test theoretical models of dynamic processes in cells.
https://doi.org/10.7554/eLife.10167.002Introduction
Recently a model was presented in which the variation of numbers of a particular type of organelle (Golgi apparatus, vacuoles or peroxisomes) observed in cells was proposed as a diagnostic indicator of the relative importance of different processes by which organelles can be formed and destroyed (Mukherji and O’Shea, 2014; see Mukherji and O’Shea, 2015 for a correction). Here we reexamine the mathematical analysis of this model and show that further insight can be gained from considering exact calculations of the equilibrium distributions. For conciseness we will refer to the model, and the analysis in the associated paper (Mukherji and O’Shea, 2014), by the abbreviation SMOP (stochastic model of organelle production).
Analysis
The SMOP model in the context of “birth and death” models
In the SMOP model, four processes are envisaged for the production and destruction of organelles: de novo synthesis, fission, fusion and decay. These four processes are characterised by one rate constant each, defined in the SMOP paper as k_{de novo}, k_{fission}, k_{fusion}, and γ. Following the definitions in the SMOP paper, the probabilities of each of the four processes occurring in the next small time period δt are given in Table 1. We also include in this table the total rate of each process that would be observed instantaneously in a large population of N cells.
Models involving processes of this type are generically termed “birth and death” processes and have a very long history of analysis in the context both of the life sciences (e.g. evolution; Yule, 1924) and in the context of physical processes (e.g. detection of cosmic rays; Furry, 1937). Accessible discussions can be found in several books (Bailey, 1990, which is a reissue of the classic text from 1964; Taylor and Karlin, 1998). In such analyses, the three processes of de novo production, production by fission, and loss by first order decay are often termed immigration, birth and death, respectively. Immigration is used for a process that increases the number of individuals but does not require any other individuals already to be present. Birth is the process by which one individual gives rise to a second individual. Death is a process by which a particular individual is lost from a population with a probability that is independent of any other members of the population. Analyses including a fusion term are much less common.
As there are a considerable number of possible combinations of the four processes that might be active, we will use a notation here to define a model by listing in curly brackets the active production processes, followed by the active destruction processes, separated by a semicolon. Any process that is not mentioned has a rate constant of zero. Thus the model with de novo, fission and decay terms would be denoted {de novo, fission; decay}.
During single cell simulations based upon the equations in Table 1 the number of organelles will fluctuate (Figure 1A) and one can ask what fraction of time, f_{n}, does a cell spend having n = 0, n = 1, n = 2, etc. organelles. This is equivalent to asking what fraction of a large ensemble of cells have n = 0, n = 1, n = 2, etc. organelles at one moment in time. In treatments of stochastic systems, the values of f_{n} would normally be described as the probability distribution for a cell having n organelles. In terms of a population of cells it can be described as a population distribution.
The simulations in Figure 1A illustrate the important point that a simulation is always started from some arbitrary starting point, and that a period of time must elapse before the simulations can be considered to be independent of this starting point. If the distribution f_{n} is evaluated at different times after the starting point of simulations (Figure 1C) then different distributions are obtained; hence the distribution is “time dependent”. As the probability (or population) distribution varies in time (and since the system can in principle be started from any state) then it is not strictly possible to talk of “the distribution” for a stochastic system. However in many birth and death models the system will settle down to a limiting distribution, independent of the starting states of the cell(s) as in the cyan and red curves in Figure 1C. Such a situation corresponds to a state of dynamic equilibrium that is familiar from chemical kinetics. Thus the terms limiting, equilibrium (or steady state) distributions can be used in this context interchangeably. The conditions for the models studied here to have limiting distributions is discussed further below; an example of a set of parameters for which a SMOP model does not yield a limiting distribution is shown in Figure 1D.
Assuming that a limiting probability distribution does exist there are two basic sampling methods by which it can be measured, irrespective of whether one is making experimental observations or performing simulations. One method is to note n at a set of time points for a single cell (such as points drawn from the trajectories in Figure 1A), and the other is to take a large number of cells at one point in time, and measure n across this ensemble of cells (Figure 1B,C). The former would require a time dependent set of observations that may be difficult to obtain experimentally. The latter approach is equivalent to making experimental observations on a large field of view of cells, or of making repeated simulations. However there is a large caveat that when the measurements are made one must be convinced that the cells have had “long enough to reach equilibrium” since the last significant perturbation to the system. If this is not the case then the distribution measured will be contaminated with contributions from nonequilibrium distributions (such as from the magenta, yellow, green and blue curves in Figure 1C). Perturbations include the choice of an arbitrary starting point in simulations, and effects such as cell division and change of growth conditions in experimental data.
Does the SMOP analysis imply equilibrium distributions?
It is implicitly assumed in the SMOP paper that the populations are to be considered to be at equilibrium (or that the probability distributions are in their limiting form) for all three cases of the analysis via simulations, from experimental data, or via the fluctuation dissipation theorem; a comment has been added to the original articles to clarify this assumption for the simulations (see the comment dated November 23, 2015 on Mukherji and O'Shea, 2014). For the experimental data this seems a reasonable assumption, although dynamic population data are really required to fully settle this issue. The fluctuation dissipation theorem method implicitly assumes a steady state (Paulsson, 2005).
Derivation of recurrence relation for the distribution of organelles in the {de novo, fission; decay, fusion} model
By applying an equilibrium condition it is straightforward to derive precise relations for the distributions in the three scenarios considered in the original SMOP paper, and hence avoid the approximations introduced by the use of the fluctuation dissipation theorem.
At equilibrium, the rate at which the population gains cells with n + 1 organelles due to cells with n organelles gaining one organelle must equal the rate at which the cells with n + 1 organelles lose one organelle. The reasoning is the same as for standard treatments of dynamic equilibrium between two states (as in a chemical reaction), and the complete justification of this when there are multiple states (i.e. cells with n = 0, n = 1, n = 2, etc., organelles) is given in Appendix 1.
Thus at equilibrium
which gives
From Equation 2 the exact distribution of organelle numbers at equilibrium can be calculated for a model involving any combination of the four processes, without recourse to random number based simulations and the attendant issues of ensuring adequate sampling precision.
An explicit numerical example of the use of Equation 2 to generate a distribution is given in Appendix 2. Briefly, an arbitrary value for f_{0} is chosen; f_{1} is then calculated from f_{0}; f_{2} is calculated from f_{1}; f_{3} is calculated from f_{2}; etc. Finally the entire distribution is normalised, which removes any dependence on the initial choice for f_{0}. Equation 2 is often termed a recurrence relation (or sometimes recursion relation or difference equation) as it allows successive terms in a distribution to be calculated from earlier terms.
Application of recurrence method to Golgi and vacuole models
The recurrence relation readily allows the derivation of precise distributions for the case of the model applied to Golgi ({de novo; decay}, Appendix 3) and vacuoles ({fission; fusion}, Appendix 4). For the Golgi, a Poisson distribution is obtained as the limiting distribution in accord with the SMOP analysis. However for vacuoles a truncated Poisson distribution is obtained, and not the shifted Poisson distribution that is reported in the SMOP analysis. Although the difference between these distributions is quite subtle (Appendix 4), the variation of Fano factor with <n> is significantly different: the Fano factor for the truncated Poisson approaches 1 much more rapidly (Figure 2, green curve) than for the shifted Poisson (Figure 2, black curve).
In Figure 2, it can be seen that the experimental values quoted in the SMOP analysis are in excellent agreement with the incorrect prediction, whilst the agreement with the corrected theoretical prediction is much less good. This greatly weakens the argument that the SMOP model makes “quantitatively accurate predictions” or that it therefore correctly accounts for the behaviour of the vacuole population.
Application of recurrence method to peroxisome models
Having established the value of analysing the SMOP model with the recurrence method, we move on to the case of peroxisomes.
Recurrence relation demonstrates that {de novo, fission; decay} yields a negative binomial distribution
For peroxisomes, the primary model discussed in the SMOP analysis is a model in which peroxisomes can potentially form both de novo and by fission, and fusion is considered to be negligible. We show in Appendix 5 that such a {de novo, fission; decay} model has a limiting distribution that is a negative binomial distribution.
As a result, exact expressions (for all parameter values) are readily obtained (Appendix 5) for the mean and Fano factor for this model,
Combining Equations 3,4 gives an alternative form for the Fano factor
This latter equation is the form given in the SMOP analysis as the k_{fusion }= 0 limit of the approximate Equation 1 of the SMOP paper. By explicitly applying the equilibrium assumption we have therefore shown that this expression is exact in this case for all values of the parameters.
Evaluation of Fano factor for the boundary marking equal de novo and fission rates
An attractive claim made for the SMOP analysis is that the value for the Fano factor can be used to distinguish between different modes of organelle production.
In Figure 3D of the original SMOP paper, a green dashed line was placed at ${\sigma}^{2}/>">n$ and stated as “marking the boundary between de novo synthesis and fission dominated organelle production”. Thus, as presented in Figure 3D of the original version of the SMOP paper, the Fano factor for cells grown on oleate (${\sigma}^{2}/>">n$ appears to lie significantly above the line where the rates are equal, and thus terms such as “fission dominated biogenesis” are used widely in the SMOP paper.
However, according to Equation 5 the correct value with equal contributions from the two processes is ${\sigma}^{2}/>">n$. A correction to this effect has now been issued. The value of ${\sigma}^{2}/>">n$ corresponds to the extreme case of zero fission and production solely by de novo production (Figure 3). In the corrected version of Figure 3D in the SMOP paper it is clear that the inferred fission rate for growth on oleate is only barely greater than the de novo synthesis rate and the term “dominance” does not seem to be appropriate. It is clear from Figure 3 that the inferred contribution from fission grows only rather slowly as the Fano factor increases above a value of two. In other words, for fission to be truly dominant a Fano factor would have to be observed that was far greater than any reported in the SMOP paper.
Restrictions on the relationships between parameters in order for limiting distributions to exist
Since it is inherent in the SMOP analysis that the distributions are limiting/equilibrium distributions, it is important to consider whether a limiting distribution will ever be reached.
In the case of the peroxisome model, {de novo, fission; decay}, for the system to have a limiting or equilibrium distribution there is a strong restriction on parameters, namely that k_{fission} < γ. To see why this restriction exists, consider first the more simple {de novo; decay} model. Considering a single cell, the number of organelles will be limited by the fact that as n grows larger then the decay process (whose rate increases proportionally to n) will become increasingly more likely than the de novo formation process (which is independent of n). If the fission process is introduced then there is now a production term that also increases proportionally to n. If k_{fission} exceeds γ then the fission process will always exceed the decay process and the number of organelles will grow without limit.
The behaviour as k_{fission} becomes similar to γ can also be seen from Equations 3,4 since if k_{fission} is increased from zero until it becomes equal to γ then <n> and ${\sigma}^{2}/>">n$ both diverge. Thus at first sight there appear to be two different ways in which the Fano factor can become very large. One way, as extensively discussed in the SMOP analysis and above, is if the mean rate k_{fission}<n> becomes large compared to the mean rate k_{de novo}. The other way, that becomes clear from our analysis, is if the rate constant k_{fission} becomes similar to the rate constant γ. The connection between these two relations is that <n> is a function of all three rate constants. This emphasises the hidden complexity of the interplay of the parameters in this model.
Discussion
Our motivation for this in depth analysis of the SMOP method was sparked by the claim that it could differentiate between fission and fusion dominated mechanisms of peroxisome biogenesis. That the distribution of numbers of organelles can give a clue to the mechanisms by which organelles are formed is a very elegant idea, and the SMOP analysis combines this idea with a very simple kinetic model. As we explored the system further we realised that the fluctuation dissipation theorem result was not necessary for analysis of the system, and that enforcing the equilibrium condition, that was implicit in the work already, greatly simplified the analysis.
There are a number of factors that cause us to question the utility of the SMOP model. A main piece of evidence for the correctness of the model was the agreement of the experimentally observed Fano factors for the vacuole data with those from the model. We have shown this agreement to be much less perfect than originally demonstrated. We have also shown that there is a strong interplay between different parameters in the model. This means that the agreement of experimental data with the model is not as compelling as originally presented and that the interpretation of experimental observations back to mechanistic conclusions is open to question. We hope that our analysis will stimulate discussion as to whether, for instance, the SMOP model captures the key features of the underlying processes and is just lacking some details; or whether the model fundamentally lacks key aspects of feedback. The assumption that observations of cells grown in batch culture faithfully report equilibrium distributions also requires further verification.
A key conclusion of the SMOP analysis is that the contribution of fission to peroxisome biogenesis is negligible (<10%) when yeast cells are grown on glucose, but "dominant" when they are grown on oleate. This is an area of some contention (Hoepfner et al, 2005; Motley and Hettema, 2007), and a recent model relied on fission of peroxisomes during organelle inheritance as the proliferation mechanism (Knoblach et al, 2013). Our analysis has shown that the term “dominant” is misleading, and that the data reported for haploid cells grown on oleate indicates approximately equal contributions from the two processes.
Nevertheless the model does suggest that the proportion of production by fission increases by about a factor of five on switching to oleate growth. Supporting evidence for this was provided by the observation of the reduction in the inferred fission contribution in cells grown on oleate in which the fission factors Vps1 or Dnm1 (or Fis1, an accessory factor of Dnm1) were deleted. On the other hand, no data were shown for glucose grown cells harbouring the same deletions. On glucose the Fano factor is reported in the SMOP analysis as 1.1, with <n> = 3, and from Equations 3,4 one obtains k_{de novo }= 2.7γ and k_{fission }= 0.1γ. The model then implies that on deletion of the fission pathway (i.e. setting k_{fission }= 0) then <n> would only drop by 10%. Peroxisome count data has been reported recently (Fig S4, Motley et al., 2015), with values of <n> = 4.9 in WT cells, <n> = 1.5 in vps1Δ cells and <n> = 1.2 in vps1Δdnm1Δ cells. The drops in peroxisome numbers in vps1Δ and vps1Δdnm1Δ cells are much greater than the 10% estimated above from the SMOP model. There is also peroxisome count data in Kuravi et al. (2006), which gives <n> = 1.6 (WT), <n> = 1.2 (vps1Δ), <n> = 1.7 (dnm1Δ), <n> = 0.9 (dnm1Δvps1Δ). The drop in <n> for the dnm1Δvps1Δ again conflicts with the idea that fission is such a small contributor to the biogenesis process. The discrepancies between these various data possibly arise from difficulties in quantifying peroxisome numbers, especially when cells contain a large number of small (and therefore low fluorescence) peroxisomes as may be the case when fission is a strong contributor to biogenesis. The problem of peroxisome counts depending on the brightness of fluorescent markers has been commented on by Jung et al. (2010).
There is a continuing push for cell biology to become more quantitative, and to be subject to the use of rigorous models as are common in the physical sciences. Such a push raises significant challenges not only in terms of developing tractable models and justifying the underlying assumptions, but also in terms of the application of the model to complex experimental data. In particular this work highlights the need for greater accounting for detection limits and intensity distributions, as well of time dependent issues, in the reporting and analysis of organelle count data, if these are to be used to infer details of organelle biogenesis mechanisms.
Materials and methods
All calculations were performed in python 2.7, running either via Cygwin under Windows 8.1, or Linux Mint 17.0. The code used for running stochastic simulations and for calculating distributions via Equation 2 is given in Source code 1,2. Time units are arbitrary. The time step for the simulations in Figure 1 was 0.0001 units. Equation 2 can be reformulated in terms of three parameters, e.g. k_{de novo}/γ, k_{fission}/γ, k_{fusion}/γ; thus for example the parameters {k_{de novo }= 2.0, k_{fission }= 0.9; γ = 1.0, k_{fusion }= 0.02} and any uniformly scaled set of parameters (e.g. {k_{de novo }= 4.0, k_{fission }= 1.8; γ = 2.0, k_{fusion }= 0.04}) yield the same limiting distribution. As in the SMOP paper, the Fano factor is defined here as $\frac{{\sigma}^{2}}{>}$.
Appendix 1: Detailed balance
Consider an ensemble of cells, where cells with for example n = 3 organelles can only arise from cells with n = 2 or n = 4 organelles.
We can represent the kinetic processes thus
R_{n=1}_{→n=0}  R_{n=2}_{→n=1}  R_{n=3}_{→n=2}  R_{n=4}_{→n=3}  R_{n=5}_{→n=4}  
n = 0  ⇆  n = 1  ⇆  n = 2  ⇆  n = 3  ⇆  n = 4  ⇆  etc. 
R_{n=0}_{→n=1}  R_{n=1}_{→n=2}  R_{n=2}_{→n=3}  R_{n=3}_{→n=4}  R_{n=4}_{→n=5} 
Where for example R_{n=2→n=3} is the rate of formation of cells with n = 3 due to cells with n = 2 gaining one organelle, and R_{n=3→n=2} is the rate of loss of cells with n = 3 due to such cells losing one organelle.
If the system is in equilibrium (i.e. the values of f_{0}, f_{1}, f_{2}, etc. are not changing) then a form of the principle of detailed balance applies. The principle of detailed balance asserts that for instance
To show this, consider first the state with n = 0. At equilibrium the rate of formation must equal the rate of loss,
Now considering the state with n = 1, the equivalent relation contains more terms
But, crucially, by substituting in from Equation 7 this simplifies to
And by proceeding along the chain the general relationship emerges that
For i = 1, 2, 3, etc.
In other words at equilibrium every individual reversible process must have balancing rates. The full “principle of detailed balance” for molecular systems contains some greater subtleties related to thermodynamics, that allow the principle to be applied to more complex networks.
Appendix 2: Recurrence example
The recurrence relation can be applied in an algebraic form as in the analysis of the various models in the paper, but can also be used to rapidly calculate numerical distributions directly.
As an example, consider the model {k_{de novo }= 2.0, k_{fission }= 0.333; γ = 1.0, k_{fusion }= 0.0333}
Equation 2 becomes
We first calculate provisional (unnormalised) values for f_{0}, f_{1}, f_{2}, etc.
First set f_{0 }= 1. This is an arbitrary choice and the final results do not depend upon it.
Using Equation 11 with n = 0,
Then using Equation 11 with n = 1, and the above value for f_{1},
And again using Equation 11 with n = 2, and the above value for f_{2},
This procedure can be repeated to calculate any number of terms.
For a model with a limiting distribution then for large n the f_{n} will become very small, for instance here f_{20 }= 0.00000008, and so it is not necessary to calculate more terms.
The calculated terms are then summed to give a normalising factor (NF),
The provisional values of f_{0}, f_{1}, f_{2}, etc. are divided by NF, which ensures that the fractional populations sum to one, giving
f_{0 }= 0.10, f_{1} = 0.20, f_{2 }= 0.23, f_{3 }= 0.19, f_{4 }= 0.13, f_{5 }= 0.08, f_{6 }= 0.04, etc.
Ideally it should be checked with a longer sequence that, to sufficient precision, the results are independent of the number of terms calculated.
Appendix 3: Application of recurrence relation to {de novo; decay} model
In the case of the Golgi, the model proposed in the SMOP analysis is a {de novo; decay} model. This model is also often termed an immigration death process. For this model Equation 2 becomes
So we have
etc., and so in general
To ensure that the sum of all f_{n} equates to one, it is simple to show from the power law definition of the exponential function that ${f}_{0}=\text{exp}({k}_{denovo}/\gamma )$, and hence
Thus the limiting distribution is a Poisson distribution as stated without proof in the SMOP paper. The mean is given by $>">n$ = k_{de novo}/γ.
Appendix 4: Application of recurrence relation to {fission; fusion} model
In the case of vacuoles, the model proposed in the SMOP analysis is a {fission, fusion} model. For this model Equation 2 becomes
Equations 16,22 are identical except for the replacement of the ratio k_{de novo}/γ by k_{fission}/k_{fusion} and therefore Equation 22 yields a form of Poisson distribution. It is not a standard Poisson distribution since, as noted in the SMOP paper, such a model cannot reach a state with n = 0 from a state with n > 0 since the fusion probability is zero when n = 1. This means that f_{0 }= 0. The rest of the distribution has the same form as for a regular Poisson distribution (Equation 23), except for a slightly different overall normalising factor. This distribution is known as a truncated (or sometimes zerotruncated) Poisson distribution (Equation 24). We have checked this analysis by direct calculation and simulation using the programs in Source code 1,2.
The relevant distributions can be defined as (Appendix 4 Figure 1)
Poisson:
Truncated Poisson:
Shifted Poisson:
Appendix 5: Application of recurrence relation to {de novo, fission; decay} model
In the case of peroxisomes, the model proposed in the SMOP analysis is a {de novo, fission; decay} model. For this model Equation 2 becomes
This has the same form as the recurrence relation for the negative binomial distribution. The negative binomial distribution is often expressed as the distribution for the number of successes in a set of trials (probability of success in each trial = p) before a certain number, r, of failures occurs, however the distribution arises in a wide variety of situations (Kotz et al, 1989; Taylor and Karlin, 1998; Scheaffer and Young, 2010; Grimmett and Stirzaker, 2001) and in general there is no restriction that r be an integer. The recurrence relation for the negative binomial distribution is,
and the negative binomial distribution has mean $>">n$ and variance ${\sigma}^{2}=\frac{pr}{{(1p)}^{2}}$.
Comparing Equations 26,27 we can equate $p={k}_{fission}/\gamma $ and $r={k}_{denovo}/{k}_{fission}$
Hence
and
Combining Equations 28,29 gives an alternative form for the Fano factor,
We also note that it is possible to derive an expression for the variance (and Fano factor) for the full {de novo, fission; decay, fusion} model as follows.
The total net rate of formation of organelles across the whole sample of cells, R_{tot}, is given by
At equilibrium R_{tot} is zero, hence we have
Making use of the relationship $>">(n>">n{)}^{2}$
we have
And thus we have the exact result that
This expression can only be applied when k_{fusion} ≠ 0, and it is not evident that it provides significant insight, however it is useful for cross checking other algebraic and computational results. We have been unable to find a useful closed expression for <n>.
References

BookThe Elements of Stochastic Processes with Applications to the Natural SciencesNew York: Wiley.

Dynamic changes in the subcellular distribution of Gpd1p in response to cell stressJournal of Biological Chemistry 285:6739–6749.https://doi.org/10.1074/jbc.M109.058552

An ERperoxisome tether exerts peroxisome population control in yeastThe EMBO Journal 32:2439–2453.https://doi.org/10.1038/emboj.2013.170

Dynaminrelated proteins Vps1p and Dnm1p control peroxisome abundance in saccharomyces cerevisiaeJournal of Cell Science 119:3994–4001.https://doi.org/10.1242/jcs.03166

Reevaluation of the role of Pex1 and dynaminrelated proteins in peroxisome membrane biogenesisThe Journal of Cell Biology 211:1041–1056.https://doi.org/10.1083/jcb.201412066

Yeast peroxisomes multiply by growth and divisionThe Journal of Cell Biology 178:399–410.https://doi.org/10.1083/jcb.200702167

Models of stochastic gene expressionPhysics of Life Reviews 2:157–175.https://doi.org/10.1016/j.plrev.2005.03.003

A mathematical theory of evolution, based on the conclusions of dr. j. c. willis, F.R.SPhilosophical Transactions of the Royal Society B: Biological Sciences 213:21–87.https://doi.org/10.1098/rstb.1925.0002
Decision letter

Vivek MalhotraReviewing Editor; The Barcelona Institute of Science and Technology, Spain
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your work entitled "A reanalysis of the stochastic model of organelle production" for peer review at eLife. Your submission has been favorably evaluated by Vivek Malhotra (Senior Editor) and two reviewers.
The reviewers have discussed the reviews with one another and the Senior Editor has drafted this decision to help you prepare a revised submission.
Specifically, the reviewers have raised a concern that your arguments are too aggressive and it would be best to tone down your statements. For example, you refer to "a number of errors" in the previous paper. The reviewers feel that this is incorrect or at least overstates the case. The reviewers agree that there was a major error in the O'Shea paper, namely the use of the shifted rather than the truncated Poisson distribution. Other than this, the mathematical results reported in the previous paper are correct as far as they go. Your manuscript shows that certain expressions in the Mukherji and O'Shea paper are exact, even though O'Shea and colleagues had reported them as approximate. This cannot be considered an error. The use of the word "dominant" is misleading, but is not a mathematical error. An abbreviated summary of the reviews follows.
A) Assessment of the manuscript:
First, I can find nothing mathematically incorrect about the present manuscript. The arguments are sound. Unlike the Mukherji and O'Shea paper, results are explicitly derived and all steps are detailed so readers are convinced of the correctness of these results. The paper is written in a rather pedagogical manner, the various results presented here are basically textbook material. For instance, deriving the detailed balance condition, explaining the stochastic simulation code, deriving the Poisson distribution, etc. would be standard graduatelevel material typically not reported in a physics paper.
I do think a good point made in the present manuscript is the issue of identifiability, i.e. is it possible to identify an underlying microscopic model based on certain macroscopic measurements? It is clear from the new analysis that many combinations of parameter values can "mix together" to produce indistinguishable Fano factors, etc. However, on its own I don't think this analysis would have been considered a strong contribution. It is only as a critique of the Mukherji and O'Shea paper that this manuscript might be considered.
B) The contribution of the Mukherji and O'Shea paper:
The Mukherji and O'Shea paper made three major contributions:
1) It put forward the idea that heterogeneity in organelle number could be used to distinguish competing models of organelle biogenesis.
2) It made predictions from different models of biogenesis.
3) It reported experiments and compared them with predictions.
C) The critique of the Mukherji and O'Shea paper:
There are two major concerns raised about the Mukherji and O'Shea paper:
1) Were the simulations and the experiments truly done so the equilibrium assumption is correct?
2) The use of the shifted rather than the truncated Poisson distribution.
There are also two minor concerns:
3) What is the basis for various approximate Fano factors reported in the Mukherji and O'Shea paper?
4) The use of the term "dominance" in many places gives the reader the wrong impression.
Let's analyse each of these in turn.
1) During the review of the Mukherji and O'Shea paper, the issue of cell division as a confounding factor was explicitly raised. The authors' response was to provide more detailed simulations in which they showed that the Fano factors did not vary greatly when partitioning of organelles during division was taken into account, and was approximately as observed when the only mode of organelle reduction was a firstorder decay term. This argument was convincing, assuming that equilibrium had truly been reached. I will add that the idea mentioned in the present manuscript (subsection “Does the SMOP analysis imply equilibrium distributions?”) that only if the time to equilibrium is short compared to the cell cycle time will the limiting distribution be valid, is a matter of scope. If the model explicitly includes cell division as one of the processes, that larger model will itself reach a limiting distribution. I agree that Mukherji and O'Shea provided no evidence that equilibrium had been reached, either in the simulations or the experiments, and were not requested to by the reviewers. This condition is typically assumed to be correct, and rarely are authors asked for evidence.
2) Is there any evidence that the simulations of the Mukherji and O'Shea paper had not reached equilibrium? I would assume that Mukherji and O'Shea would have done the correct validation steps, and there is some evidence of that in their Figure 1 for example. However, it is then difficult to understand how the error about the truncated vs. the shifted Poisson distribution could have arisen. I admit that this point could have been discovered by the reviewers if they had explicitly repeated some of the reported calculations in the Mukherji and O'Shea paper, and I for one overlooked it. It is unclear how the exact stochastic simulations could be in agreement with an incorrect shifted Poisson prediction. In the Mukherji and O'Shea paper no simulation matching the shifted Poisson prediction is given. Perhaps this simulation was never done, otherwise Mukherji and O'Shea would have discovered the discrepancy.
3) For certain expressions in the Mukherji and O'Shea paper, which arise out of an approximate calculation using the fluctuationdissipation theorem (subsection “Incorrect analysis of cases with k_{fusion} ≠ 0” of this manuscript), it is not explained how they were derived. This is a valid but minor concern; Mukherji and O'Shea would likely be able to provide the derivation. The point that the approximations were not necessary given the existence of an exact result is worth reporting but does not on its own detract from the original findings. In particular, many expressions derived in the Mukherji and O'Shea paper are correct, as shown in the present manuscript.
4) The misleading use of the word "dominant". There are scientific contexts where "dominance" is interchangeable with "greater than". E.g. in game theory, the word "dominance" is applied when the payoff from one strategy is just greater than from other strategy. Mukherji and O'Shea used the word "dominance" or "dominates" throughout their paper, but "contributed greater than" or some such term might have been more appropriate. I agree that this gives the reader the wrong impression, but does not, on its own, invalidate the results of Mukherji and O'Shea.
D) Is the critique valid?
I had said in my original review and I reiterate there: the ideas in the Mukherji and O'Shea paper were worth reporting, and I had hoped the original Mukherji and O'Shea paper would stimulate discussion and be improved upon with new analyses and measurements. The present manuscript in a sense validates this hope: discussion has indeed been stimulated. However, it is true that the original paper contained one major error not caught in the review process: use of the shifted rather than the truncated Poisson distribution. The present author argues that this weakens the central attraction of the previous paper, which was the strong agreement of theory and experiment. I would argue that the agreement itself, which was admittedly a striking feature, was not the central point on which the Mukherji and O'Shea paper stood. If the original paper had reported the correct distribution and highlighted the discrepancy with measured data, that would have also been food for thought.
Overall, I feel this is a more accessible and comprehensive theoretical analysis than the Mukherji and O'Shea paper provided, and serves to continue the original discussion.
https://doi.org/10.7554/eLife.10167.010Author response
Specifically, the reviewers have raised a concern that your arguments are too aggressive and it would be best to tone down your statements.
I have softened the language throughout the text. In addition, some of the references to parts of Figure 1 and Appendix 4 Figure 1 have also been improved.
I have also taken the opportunity to add one sentence where I explain how I hope this will stimulate discussion of the underlying premises of the model (Discussion: “Our analysis has shown that the term […] contributions from the two processes”).
For example, you refer to "a number of errors" in the previous paper. The reviewers feel that this is incorrect or at least overstates the case. The reviewers agree that there was a major error in the O'Shea paper, namely the use of the shifted rather than the truncated Poisson distribution. Other than this, the mathematical results reported in the previous paper are correct as far as they go.
I have toned down some of these statements.
However the fact that the value of the Fano factor given by the expression 1 k_{fission}<n>/k_{de novo} is equal to 2 when the terms k_{fission} <n> and k_{de novo} are equal and not to 1 as stated in the SMOP paper is a clear second mathematical error and was the reason why the k_{fission} term appeared so “dominant” in Figure 3D, since the value for oleate growth is so close to 2 (almost within the stated error).
Your manuscript shows that certain expressions in the Mukherji and O'Shea paper are exact, even though O'Shea and colleagues had reported them as approximate. This cannot be considered an error. The use of the word "dominant" is misleading, but is not a mathematical error. The detailed comments of the reviewers are pasted below.
I have used the term “misleading” now so that I hope that will lead people to look more closely at the original arguments – and that that can lead to some of the healthy discussion that the reviewers refer to.
I am particularly keen to know the origin of the expression
$\frac{{\sigma}^{2}}{>}$since I cannot find a way that it follows from their Equation 1 but maybe they have derived it via a different route. I am very happy to remove that comment if necessary, as the peroxisome model plus fusion does not seem to have much relevance given prevalent models of peroxisome biogenesis (detailed query is below).
A) Assessment of the manuscript:
First, I can find nothing mathematically incorrect about the present manuscript. The arguments are sound. Unlike the Mukherji and O'Shea paper, results are explicitly derived and all steps are detailed so readers are convinced of the correctness of these results. The paper is written in a rather pedagogical manner, the various results presented here are basically textbook material. For instance, deriving the detailed balance condition, explaining the stochastic simulation code, deriving the Poisson distribution, etc. would be standard graduatelevel material typically not reported in a physics paper.
This area is so multidisciplinary that I felt I should make the Discussion as accessible to as many workers in the field as possible. In addition, giving explicit examples help remove as much ambiguity as possible from definitions. I give the derivation of the Poisson distribution to show that it does arise from the recurrence relation: many people will be much more familiar with the Poisson distribution in the context of the number of successes of a large number of lowprobabilityofsuccess trials, which is really rather different from the birthdeath type situation. I also give the Poisson (Golgi) case to show that my analysis does agree with the SMOP analysis in that case.
I give a full explanation of the principle of detailed balance as used here because the application here seems to me to be different from the thermodynamic principle of detailed balance. In a cellular system it seems there is no objection to a cycle having stable populations and yet sustaining a net flux round the cycle, whereas that is forbidden in closed thermodynamic systems. In the case of the present model it is the linear (and n>=0) nature of the set of linked reversible processes that lead to detailed balance. If there are accessible presentations of this material readily available and couched in sufficiently identical terms then I am more than happy to direct readers to those instead.
I do think a good point made in the present manuscript is the issue of identifiability, i.e. is it possible to identify an underlying microscopic model based on certain macroscopic measurements? It is clear from the new analysis that many combinations of parameter values can "mix together" to produce indistinguishable Fano factors, etc. However, on its own I don't think this analysis would have been considered a strong contribution. It is only as a critique of the MO paper that this manuscript might be considered.
In light of this comment, and the encouragement that this paper can form part of a healthy discussion, I have added a short paragraph discussing the implications of the underlying assumptions of the model (Discussion: “Our analysis has shown that the term […] contributions from the two processes”).
[…]
C) The critique of the Mukherji and O'Shea paper:There are two major concerns raised about the Mukherji and O'Shea paper: 1) Were the simulations and the experiments truly done so the equilibrium assumption is correct? 2) The use of the shifted rather than the truncated Poisson distribution. There are also two minor concerns: 3) What is the basis for various approximate Fano factors reported in the Mukherji and O'Shea paper? 4) The use of the term "dominance" in many places gives the reader the wrong impression.
I feel this list omits the crucial miscalculation of the value of the Fano factor for the case where the fission and de novo rates are equal in the peroxisome model, as I have discussed above.
Let's analyse each of these in turn.
1) During the review of the Mukherji and O'Shea paper, the issue of cell division as a confounding factor was explicitly raised. The authors' response was to provide more detailed simulations in which they showed that the Fano factors did not vary greatly when partitioning of organelles during division was taken into account, and was approximately as observed when the only mode of organelle reduction was a firstorder decay term. This argument was convincing, assuming that equilibrium had truly been reached. I will add that the idea mentioned in the present manuscript (subsection “Does the SMOP analysis imply equilibrium distributions?”) that only if the time to equilibrium is short compared to the cell cycle time will the limiting distribution be valid, is a matter of scope.If the model explicitly includes cell division as one of the processes, that larger model will itself reach a limiting distribution. I agree that Mukherji and O'Shea provided no evidence that equilibrium had been reached, either in the simulations or the experiments, and were not requested to by the reviewers. This condition is typically assumed to be correct, and rarely are authors asked for evidence.
The cell division model is only included as a brief supplement to Figure 1 of the SMOP paper, and only for the Golgi ({de novo; fission}) model. Almost no details are given and that is significant because the issue of timescale becomes relevant in such a model. How does the “regular time interval” relate to the timescale of the various processes?
2) Is there any evidence that the simulations of the Mukherji and O'Shea paper had not reached equilibrium? I would assume that Mukherji and O'Shea would have done the correct validation steps, and there is some evidence of that in their Figure 1 for example. However, it is then difficult to understand how the error about the truncated vs. the shifted Poisson distribution could have arisen. I admit that this point could have been discovered by the reviewers if they had explicitly repeated some of the reported calculations in the Mukherji and O'Shea paper, and I for one overlooked it. It is unclear how the exact stochastic simulations could be in agreement with an incorrect shifted Poisson prediction. In the Mukherji and O'Shea paper no simulation matching the shifted Poisson prediction is given. Perhaps this simulation was never done, otherwise Mukherji and O'Shea would have discovered the discrepancy.
3) For certain expressions in the Mukherji and O'Shea paper, which arise out of an approximate calculation using the fluctuationdissipation theorem (subsection “Incorrect analysis of cases with k_{fusion} ≠ 0” of this manuscript), it is not explained how they were derived.This is a valid but minor concern; Mukherji and O'Shea would likely be able to provide the derivation.
If it can be obtained I am very happy to drop this point. However I cannot see that it can follow from Equation 1. Equation 1 is
$\frac{{\sigma}^{2}}{>}$If the limit is looked for under the condition k_{fusion}<(n(n1))>) ~ γ <n>, it is not clear how to proceed because one immediately has terms in <n^{2}>. (It appears that one can make some progress by then using variance=<n^{2}><n>^{2} but that does not seem to lead to the quoted result).
If one instead writes the first fraction in the denominator of the left hand side of the above equation as X and try to solve simultaneously,
$\frac{{\sigma}^{2}}{>}$along with the quoted result
$\frac{{\sigma}^{2}}{>}$then one obtains an expression for X that depends on k_{fission}<n> and k_{de novo} so again this route does not seem to lead to anywhere. However as <n> is inextricably linked to all parameters then maybe there is a hidden relationship that leads to the quoted result.
The point that the approximations were not necessary given the existence of an exact result is worth reporting but does not on its own detract from the original findings. In particular, many expressions derived in the Mukherji and O'Shea paper are correct, as shown in the present manuscript.
I agree, but if quantitative methods are to be adopted more within the mainstream then I think it behoves the modelers to make sure they keep things as simple as possible, otherwise it becomes the preserve of applied mathematician “wizards”. The recurrence method is much simpler than the FDT approach and should be accessible to more of the community.
4) The misleading use of the word "dominant". There are scientific contexts where "dominance" is interchangeable with "greater than". E.g. in game theory, the word "dominance" is applied when the payoff from one strategy is just greater than from other strategy. Mukherji and O'Shea used the word "dominance" or "dominates" throughout their paper, but "contributed greater than" or some such term might have been more appropriate. I agree that this gives the reader the wrong impression, but does not, on its own, invalidate the results of Mukherji and O'Shea.
I wrestled with this and sought the opinions of colleagues about what they would understand by “dominant”. If one rate were dominant over another then some people responded by saying it would imply a condition as extreme as differing by an order of magnitude. Personally I would want to see a factor of about three. Since the SMOP paper is a quantitative paper one would expect to see a statement of the ratio of the two rates, rather than the loose word dominant. In the present context the argument is made more complex by the use of a Fano factor of 1 rather than 2 as the boundary between equal rates. This is a significant difference when trying to evaluate the meaning of a Fano factor of 2.4 /0.2. According to my analysis the rates are in the ratio 0.58:0.42, or roughly 1.4:1. Given the experimental and theoretical uncertainties this does not look to me like a secure case of dominance.
I agree that there are cases where dominant can mean “one just greater than the other”, and the game theory use is such a case. In politics an election vote distributed 0.58:0.42 would tend to be seen as indicating the dominance of a party.
Indeed if two rates act against each other then only a tiny “dominance” can lead to run away behaviour.
But here we are evaluating the relative magnitude of two rates that act in the “same direction” so I think that the requirement that dominance indicates a clear cut difference is valid.
I have amended my wording to “misleading” as suggested.
D) Is the critique valid?
A4F1I had said in my original review and I reiterate there: the ideas in the Mukherji and O'Shea paper were worth reporting, and I had hoped the original Mukherji and O'Shea paper would stimulate discussion and be improved upon with new analyses and measurements. The present manuscript in a sense validates this hope: discussion has indeed been stimulated. However, it is true that the original paper contained one major error not caught in the review process: use of the shifted rather than the truncated Poisson distribution. The present author argues that this weakens the central attraction of the previous paper, which was the strong agreement of theory and experiment. I would argue that the agreement itself, which was admittedly a striking feature, was not the central point on which the Mukherji and O'Shea paper stood. If the original paper had reported the correct distribution and highlighted the discrepancy with measured data, that would have also been food for thought.
The extremely close agreement of the observed and (incorrectly) predicted Fano factors seems to have provided a key link in the argument that the underlying model is a good one. That, in itself, is part of a fairly standard inductive scientific method with its attendant strengths and weaknesses. There was very little justification otherwise of the underlying model, and certain aspects of the model deserve more discussion.
For instance the use of a rate for fusion proportional to n(n1) would seem appropriate if the organelles could be considered to be freely moving and well mixed as the molecules within a liquid. But whether it is appropriate for objects such as vacuoles seems harder to justify. It implies that the likelihood (within the next small deltat) of a vacuole fusing with another when n = 3 is double that when n = 2. One could envisage that when n = 2 there could be a very similar degree of total contact area with other vacuoles as when n = 3, and so the fusion probability would be very similar. It would depend on vacuole size, composition and “packing”. I am not aware that sufficiently detailed analysis of such factors exists: it is certainly not quoted in the SMOP paper.
The observation of very good agreement with the model would reduce the impact of such an objection. However, in the absence of such agreement, then, such discussion is vital: is the model basically a good one but just lacking a little detail or is the model fundamentally lacking key aspects of feedback? I hope this will stimulate more discussion on these issues.
https://doi.org/10.7554/eLife.10167.011Article and author information
Author details
Funding
No external funding was received for this work.
Acknowledgements
I thank Ewald Hettema, Paul Galvin and Ian Sudbery for helpful discussions, and Eleanor Stillman for pointing me towards helpful textbooks on stochastic processes.
Reviewing Editor
 Vivek Malhotra, The Barcelona Institute of Science and Technology, Spain
Publication history
 Received: July 17, 2015
 Accepted: August 26, 2015
 Version of Record published: January 14, 2016 (version 1)
Copyright
© 2016, Craven et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,023
 Page views

 131
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.