In a seminal article in virology ('Sur l’unité lytique du bactériophage' Comptes rendus des séances de la Societé de biologie et de ses filiales, 1939, 130, pp. 904-907) the Nobel prize winner Salvador Luria applied the Poisson distribution to demonstrate that one bacteriophage particle was sufficient to lyse a culture of 80 million bacteria.

He defined the probability P that the actual number of virions in a tube was x as:

P = (nˣ/x!)e⁻ⁿ

where n is the average number of virions in a series of tubes. He then defined the probability Q that x would be less than m as the sum between x=0 and m-1 of P. For m=1:

Q = e⁻ⁿ

He then provided a set of numbers of virions x and the corresponding set of Q probabilities experimentally determined:

N = c(1, 2, 4, 8, 15)
Q = c(0.736, 0.586, 0.357, 0.135, 0.022)

However, I do not understand what actually goes into the P equation and how to derive Q. I have defined a Poisson function, considering the number of particles x as a success, but I do not understand where to get the number of trials (Luria used 100 tubes, but pois_f(X, 100) does not work either). I used a function to get a vector of probabilities:

pois_f = function(x, n) {        
  p = ((n^x)/factorial(x)) * exp(-n)
p1 = pois_f(N, length(N))

I tried to get a Q by exponentiating the negative number of particles found, but the slope is off:

p2 = exp(-N)

enter image description here

If p1 is wrong because I am measuring the poison distribution (assuming that I have used the function correctly, perhaps I need to invert the terms...) instead of Q; p2 should be Q, but it is also off the target.

How do I correctly calculate Q? What measurements are required to define Q? Is Q 1-P?

PS: for reference, here as some extracts of the original paper: Definition of P (poison distribution): enter image description here

Definition of Q (cumulative distribution): enter image description here

Original graph (data here based on lower set, I): enter image description here

New graph based on the updated information: enter image description here


1 Answer 1


Your P equation is just a probability mass function of Poisson distribution, but what you need is the cumulative distribution function for the distribution. They are both available in R as dpois and ppois respectively.

What you are missing is the $\lambda$ parameter (N in your case). It is not the number of trials, but rather something like the average number of virions (I didn't read the paper, so take it with a grain of salt). It definitely is not the number of tubes. To replicate the result you would need to know the parameter (either know what it was or how it was calculated by the author).

Your comment has shed some light on what you are describing. To reproduce the result, you need to show the cumulative distribution function of the Poisson distribution parametrized by $\lambda = N$. The cumulative distribution function $F$ calculates the probability that the random variable $X$ is less or equal to some value $x$. In your case, you are interested in the "probability $Q$ that $x$ would be less than $m$" since the Poisson distribution is discrete and non-negative, it's just $P(X < 1) = P(X = 0)$. In the case of Poisson distribution,

$$ P(X = x) = \frac{\lambda^x e^{-\lambda}}{x!} $$

for $x$ equal to $0$ reduces to $e^{-\lambda}$, as in the paper. This means that what they are calculating is $P(X=0) = e^{-\lambda}$ for different values of $\lambda$, presumably equal to N.

This is consistent with what you described, but I have no idea why the numbers are not consistent with the paper. I don't have access to the paper and don't know French, so cannot comment beyond what was described by you.

  • $\begingroup$ Thank you for your reply. The article is in French and I don't know how to attach it here as a reference. Anyway, on second reading, what I called X array is actually N: the expected number of virions as opposed to X, the actual number of virions. I have updated the post accordingly. The cumulative distribution function is Q, or better: Qₘ(n), which is the sum from x=0 to x=m-1 of P. For m=1, Qₘ(n) = e⁻ⁿ ≡ exp(-λ), which leaves us with the computation of p2. Essentially, λ≡N. Yet, the slope is off... $\endgroup$
    – Gigiux
    Commented Apr 3, 2023 at 10:21
  • $\begingroup$ @Gigiux thanks for comment and see my edit. $\endgroup$
    – Tim
    Commented Apr 3, 2023 at 10:39
  • $\begingroup$ Thank you. The problem with ppois is that it introduces an additional function not present in the original paper. I could not paste the all paper but I extracted the formulae I used for the computation. There was also a typo in the definition of pois_f which I have corrected. Yet, the simulation even with ppois(1, N) fit poorly, especially in log space... $\endgroup$
    – Gigiux
    Commented Apr 4, 2023 at 12:35
  • $\begingroup$ @Gigiux "probability Q that x would be less than m" is the ppois function, there is no other cumulative distribution function for Poisson distribution. If you want to implement it manually, feel free to do so, but it's already there. I don't know what you mean by getting incorrect results from the simulations. $\endgroup$
    – Tim
    Commented Apr 4, 2023 at 12:45
  • $\begingroup$ So I should get the same results. The 'incorrect results' is shown in the final figure: the fit is poor at high N $\endgroup$
    – Gigiux
    Commented Apr 4, 2023 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.