Counting single cells and computing their heterogeneity: from phenotypic frequencies to mean value of a quantitative biomarker

Abstract

This tutorial presents a mathematical theory that relates the probability of sample frequencies, of M phenotypes in an isogenic population of N cells, to the probability distribution of the sample mean of a quantitative biomarker, when the N is very large. An analogue to the statistical mechanics of canonical ensemble is discussed.

Keywords

large deviation principle / chemical kinetics / Boltzmann’s law / variational Bayesian method / maximum entropy principle

INTRODUCTION

Statistical analyses of data and stochastic models of mechanisms are two very different, but complementary approaches in biological research. While the former obtains a quantitative representation of high-throughput measurements [¹], the latter can provide “laws of nature” through limit theorems [²], widely called emergent phenomenon. A case in point is the theory of phase transition [³] which shows that a nonlinear stochastic dynamical system with bistability and cusp catastrophe, in the limit of time t\to \infty followed by system’s size V\to \infty, necessarily exhibits a discontinuous transition [⁴]. Another example is the recent work [⁵] which demonstrates that Gibbsian equilibrium chemical thermodynamics can be reformulated as a limit theorem in a mesoscopic chemical kinetic system, with N species and M reversible stochastic elementary reactions, as the system’s size becoming macroscopic.

With the rise of single-cell biology, one naturally is interested in the limiting behavior of the phenotypic frequencies among a population of cells, usually based on one, or several biomarkers. In this case, there is actually a very powerful mathematical result that is widely known to probabilists and statistical physicists. In this tutorial, we give an introduction of this theory and discuss its broader implications.

CHARACTERIZING HETEROGENEITY IN SINGLE CELLS

Asymptotic probaility distribution for sample frequencies of cellular phenotypes

To study phenotypic heterogeneity, let a population of N isogenic cells as independent and identically distributed (i.i.d.) realizations of random events from a set �=\{1,2,...,M\}: There are totally M possible phenotypes. Among the N cells, let n_k^{}denotes the random number of cells in the k^{\mathrm{th}} state: n_1+n_2+···+n_M=N. By phenotypic frequency, we mean f_k^{(N)}\equiv {\displaystyle \frac{n_k}{N}}.

Let p_k^{} denote the probability of a cell in the k^{\mathrm{th}} state. Then the probability distribution for the observed frequency \overrightarrow{f}=(f_1,...,f_M) being \mathrm{x}=(x_1,...,x_M) follows a multinomial distribution

(1)\mathrm{Pr}\{{\overrightarrow{f}}^{(N)}=\mathrm{x}\}={\displaystyle \frac{N!}{(N{x}_{\mathrm{1}})!(N{x}_{\mathrm{2}})!...(N{x}_M)!}}{p}_{\mathrm{1}}^{N{x}_{\mathrm{1}}}{p}_{\mathrm{2}}^{N{x}_{\mathrm{2}}}...p_M^{N{x}_M}.

Since usually N is very large in a high-throughput single-cell experiment, one can safely approximate Eq. (1) using Stirling’s formula and obtain:

(2)\begin{array}{l}\ln \hspace{0.17em}\mathrm{Pr}\{f^{(N)}=\mathrm{x}\}\\ \\ \\ =\ln \hspace{0.17em}\left({\displaystyle \frac{N!}{(N{x}_{\mathrm{1}})!(N{x}_{\mathrm{2}})!\cdots (N{x}_M)!}}{p}_{\mathrm{1}}^{N{x}_{\mathrm{1}}}{p}_{\mathrm{2}}^{N{x}_{\mathrm{2}}}...p_M^{N{x}_M}\right)\\ \\ \\ \simeq -N\left[{\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k\ln \left({\displaystyle \frac{x_k}{p_k}}\right)\right].\end{array}

Therefore, one has the asymptotic limit

(3)\phi (\mathrm{x})=-\underset{N\to \infty }{\lim }{\displaystyle \frac{\mathrm{1}}{N}}\ln \hspace{0.17em}\mathrm{Pr}\left\{{\overrightarrow{f}}^{(N)}=\mathrm{x}\right\}={\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k\ln \left({\displaystyle \frac{x_k}{p_k}}\right).

In the theory of large deviations of probability, this is known as Sanov’s theorem [⁶]. Since \phi (\mathrm{x})>\mathrm{0} except when \mathrm{x}=(p_1,...,p_M), in the limit of N\to \infty, the probability of f_k^{(N)}\ne p_k is zero, and the probability of f_k^{(N)}=p_k is one. The frequency yields the probability for an infinitely large number of i.i.d. samples. Furthermore, Eq. (2) shows that (p_1,p_2,...,p_M) are the most probable sample frequencies for a finite but large N.

Asymptotic distribution for the mean value of a biomarker

Eqs. (1) and (2) give the probability for the frequencies within the N cells distributed among the M phenotypic states. We now consider a specific biomarker \mathit{g}, which is assumed to be a well defined real-valued function of the phenotype of a cell: \mathit{g}=g_k when a cell is in the k^{\mathrm{th}} state.

It is very clear that if one knows the frequencies {\overrightarrow{f}}^{(N)}, then the mean value for \mathit{g} over the entire population of the N cells is determined:

(4){\overline{\mathit{g}}}^{(N)}={\displaystyle \frac{n_{\mathrm{1}}{g}_{\mathrm{1}}+n_{\mathrm{2}}{g}_{\mathrm{2}}+\cdots +n_M{g}_M}{N}}={\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{f}_k^{(N)}{g}_k;

since the frequencies f_k^{(N)} are random, so is {\overline{\mathit{g}}}^{(N)}. Then when N\to \infty, one expects the {\overline{\mathit{g}}}^{(N)} approaching to the expected value E[\mathit{g}]. This is easy to show:

(5)\underset{N\to \infty }{\lim }{\overline{\mathit{g}}}^{(N)}={\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}\left(\underset{N\to \infty }{\lim }{f}_k^{(N)}\right)g_k={\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^m}{p}_k{g}_k=E[\mathit{g}].

What is the probability distribution for {\overline{\mathit{g}}}^{(N)} when N is very large but not infinite? One can calculate this:

(6)\begin{array}{l}\mathrm{Pr}\{{\mathit{g}}^{(N)}=y\}={\displaystyle \mathrm{\sum }_{\{\mathrm{x}:{\mathrm{\sum }}_{k=\mathrm{1}}^M{x}_k{g}_k=y\}}}\mathrm{Pr}\{f^{(N)}=\mathrm{x}\}\\ \\ \\ \simeq \exp \left(-N\underset{\{\mathrm{x}:{\mathrm{\sum }}_{k=\mathrm{1}}^M{x}_k{g}_k=y\}}{\inf }\phi (\mathrm{x})\right)\\ \\ \\ \mathrm{as}\hspace{0.17em}N\to \infty .\end{array}

We obtain Eq. (6) because among the many sets of x that give the same value y, each has a probability of e^{-N_{\phi (x)}}. Therefore, as N\to \infty, only the set with the smallest \phi (\mathrm{x}) matters. Eq. (6) indicates that for very large N , the probability distribution for the mean value of the biomarker {\overline{\mathit{g}}}^{(N)} has the form e^{-N_{\psi (y)}}, in which

(7)\psi (y)\equiv -\underset{N\to \infty }{\lim }{\displaystyle \frac{\mathrm{1}}{N}}\ln \hspace{0.17em}\mathrm{Pr}\{{\overline{\mathit{g}}}^{(N)}=y\}=\underset{\{\mathrm{x}:{\mathrm{\sum }}_{k=\mathrm{1}}^M{x}_k{g}_k=y\}}{\inf }\phi (\mathrm{x}).

In the theory of large deviations of probability, this result is known as contraction principle [⁶]. \psi (y) and \phi (\mathrm{x})are called a level-1 and a level-2 large deviations rate functions, respectively.

From phenotypic frequencies to biomarker mean values

The right-hand-side of Eq. (7) can be further carried out; this is a problem of constrained minimization using multivariate calculus:

(8)\{\begin{array}{l}\underset{\mathrm{x}}{\min }\left\{\left({\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k\ln {\displaystyle \frac{x_k}{p_k}}\right)\right\},\\ \\ {\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k{g}_k=y,\\ \\ {\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k=\mathrm{1.}\end{array}

Introducing Langrage multipliers for Eq. (8),

L(\mathrm{x},\beta ,\phi )=\left({\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k\ln {\displaystyle \frac{x_k}{p_k}}\right)+\beta \left({\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k{g}_k-y\right)

(9)+\lambda \left({\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k-\mathrm{1}\right).

Then we can find {\mathrm{x}}^{\ast }=(x_{\mathrm{1}}^{*},x_2^{*},...,x_M^{*}), {\beta }^{\ast }, and {\lambda }^{\ast } as the solution of

(10){\displaystyle \frac{\partial L(\mathrm{x},\beta ,\lambda )}{\partial x_j}}={\displaystyle \frac{\partial L(\mathrm{x},\beta ,\lambda )}{\partial \beta }}={\displaystyle \frac{\partial L(\mathrm{x},\beta ,\lambda )}{\partial \lambda }}=\mathrm{0.}

That is,

(11a)x_k^{*}={\displaystyle \frac{p_k{e}^{-{\beta }^{*}{g}_k}}{{\displaystyle {\mathrm{\sum }}_{j=\mathrm{1}}^M}{p}_j{e}^{-{\beta }^{*}{g}_j}}},

(11b)y={\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^M}{x}_j^{*}{g}_j=-{\left[{\displaystyle \frac{\partial }{\partial \beta }}\ln {\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^M}{p}_j{e}^{-{\beta }_{g_j}}\right]}_{\beta ={\beta }^{*}},

in which {\beta }^{\ast } is a function of y through Eq. (11b), which gives the function implicitly. We therefore obtain

\begin{array}{l}\psi (y)=\phi ({\mathrm{x}}^{\ast })={\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k^{*}\ln \left({\displaystyle \frac{e^{-{\beta }^{*}{g}_k}}{{\displaystyle {\mathrm{\sum }}_{j=\mathrm{1}}^M}{p}_j{e}^{-{\beta }^{*}{g}_j}}}\right)\\ \\ \\ =-{\beta }^{\ast }{\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{x}_k^{*}{g}_k-\ln {\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^M}{p}_j{e}^{-{\beta }^{*}{g}_j}\end{array}

(12)=-{\beta }^{\ast }y+{\beta }^{\ast }F({\beta }^{\ast }),

where

(13)F(\beta )=-{\displaystyle \frac{\mathrm{1}}{\beta }}\ln \hspace{0.17em}Z(\beta ),\hspace{0.17em}Z(\beta )\equiv {\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^M}{p}_j{e}^{-{\beta }^{*}{g}_j},

and {\beta }^{\ast }(y) solves \mathrm{d}[\beta F(\beta )]/\mathrm{d}\beta =y.

The above computation tells us that if one knows the values of a biomarker for all the M states of a cell, g_1,g_2,...,g_M, together with a prior knowledge of p_1,p_2,...,p_M, one should construct the Z(\beta ) function and calculate the F(\beta ) given in Eq. (13). Then the probability distribution for the mean value of the biomarker is going to be:

(14)\mathrm{Pr}\{{\overline{\mathit{g}}}^{(N)}=y\}\propto {[e^{-N\beta \{F(\beta )-y\}}]}_{\beta ={\beta }^{\ast }(y)}.

It also tells us that if one observes the mean biomarker value being \widehat{y}, then the most probable phenotypic frequencies will have a posterior form that deviates from its prior \{p_k\}:

(15)f_k\{\mathrm{conditioned}\hspace{0.17em}\mathrm{on}\hspace{0.17em}\overline{\mathit{g}}=\widehat{y}\}={\left[{\displaystyle \frac{p_k{e}^{-\beta g_k}}{Z(\beta )}}\right]}_{\beta ={\beta }^{\ast }(\widehat{y})}.

Both Eqs. (14) and (15) suggest that the functional relationship y=\mathrm{d}[\beta F(\beta )]/\mathrm{d}\beta, between the mean value of the biomarker y=\overline{\mathit{g}} and the Lagrangian multiplier \beta, or its inverse form \beta ={\beta }^{\ast }(y), are very fundamental to the probabilistic problem, in the limit of infinite sample size N\to \infty.

BEYOND AN i.i.d. POPULATION

We derived the expression in Eq. (3) based on the assumption of a population of N cells that are i.i.d. samples of a single M-state random individual with probability \{p_k\}. When there are cell-cell interactions among the individuals within a population, the mathematics immediately becomes much more involved.

Two types of research go beyond an i.i.d. population in the stochastic modeling; they were originally motivated, respectively, by chemical kinetics in solution [⁷] and Ising model for ferromagnetism of solid [⁸,⁹]. In chemical kinetics, rapid spatial movement of all “individual molecules” in an aqueous solution leads to the assumption that every individual collides with every other invidivual, and certain “reactions” can occur randomly. The Gibbs function in chemical thermodynamics is precisely like the \phi (\mathrm{x}) function in Eq. (3), for complex chemical reaction systems in equilibrium [⁵]. Actually there is a general equation, first discovered by [¹⁰], whose solution can provide \phi (\mathrm{x}) for non-i.i.d. populations. For J=M-1 reversible unimolecular reactions among M species, X_1,X_2,...,X_M, with concentrations \mathrm{x}=(x_1,x_2,...,x_M) and arbitrary non-negative functions R_{\pm j}(x) being the rates of the j^{\mathrm{t}\mathrm{h}} reaction between species j and the species M ,

(16)X_{\mathrm{1}}\underset{R_{+\mathrm{1}}}{\overset{R_{-\mathrm{1}}}{\rightleftharpoons }}{X}_M,X_{\mathrm{2}}\underset{R_{+\mathrm{2}}}{\overset{R_{-\mathrm{2}}}{\rightleftharpoons }}{X}_M,...,X_{M-\mathrm{1}}\underset{R_{+J}}{\overset{R_{-J}}{\rightleftharpoons }}{X}_M,

the equation reads

(17){\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^J}{R}_{+j}(\mathrm{x})[\mathrm{1}-e^{\partial \phi /\partial x_j-\partial \phi /\partial x_M}]+R_{-j}(\mathrm{x})[\mathrm{1}-e^{-\partial \phi /\partial x_j+\partial \phi /\partial x_M}]=\mathrm{0.}

If R_{+j}(\mathrm{x})=q_j{x}_M and R_{-j}(\mathrm{x})= r_j{x}_j, then the solution to Eq. 17 recovers the Eq. 3,

(18)\phi (\mathrm{x})={\displaystyle \mathrm{\sum }_{m=\mathrm{1}}^M}{x}_m\ln \left({\displaystyle \frac{x_m}{p_m}}\right),

in which the p’s are functions of q’s and r’s,

(19)p_m={\displaystyle \frac{{\displaystyle \frac{q_m}{r_m}}}{{\displaystyle \frac{q_{\mathrm{1}}}{r_{\mathrm{1}}}}+\cdots +{\displaystyle \frac{q_{M-\mathrm{1}}}{r_{M-\mathrm{1}}}}+\mathrm{1}}}.

The particular set of R_{\pm j}(x) represents chemical reactions in an ideal solution. A reader who had a course on freshman chemistry might recognize Eq. (18) as

G(\mathrm{x})={\displaystyle \mathrm{\sum }_{m=\mathrm{1}}^M}{x}_m{\mu }_m,{\mu }_m(x_m)={\mu }_m^o+RT\ln x_m,

where {\mu }_m is the chemical potential of m^{\mathrm{th}} specie with mole fraction x_m(not molarity) in an ideal solution, and {\mu }_m^o=-RT\ln \hspace{0.17em}{p}_m. Then \mathrm{\Delta }{\mu }_{ij}^o={\mu }_i^o-{\mu }_j^o=-RT\ln \hspace{0.17em}(p_i/p_j), where (p_i/p_j) is the equilibrium constant between species i and j [¹¹]. Apart from the RT , the Gibbs energy function is a consequence of statistical counting, which has very little to do with the energy of the atoms in the molecules [¹²].

In the second type, Ising model and alike, “individual atoms” are located at fixed lattice points in a solid, each one only interacts with its neighbours. The limit of N\to \inftyof such an interacting particle system is known as hydrodynamic limit of the stochastic model.

Cell-cell interactions in a tissue or in a culture medium can have both types: When an interaction is mediated by rapidly diffusing small molecular factors, one can safely assume the interaction is between every two individual cells in a population. If an interaction between nearby cells is mediated by slowly diffusing molecules, or due to direct contacts via mechanical interactions, gap junctions, or synapses, then a lattice model is more appropriate. Combining these two types of mathematical descriptions leads to the “reaction diffusion” paradigm [¹³] which serves the foundation for describing living phenomena [¹⁴].

DISCUSSION

Statistical mechanics and Boltzmann’s law

A reader who had a course on statistical mechanics [¹⁵] will certainly recognize Z(\beta ), F(\beta ), and {\beta }^{-1}in Eq. (13) as partition function, Helmholtz free energy, and temperature, if one identifies g_kas the energy of the k^{\mathrm{t}\mathrm{h}} state of a mechanical system. Eq. (12) then shows that F(\beta )= y+{\beta }^{-1}\psi (y)where -\psi (y) should be identified as “entropy” of the mechanical system with energy y; and it is related to F(\beta ) through a Legendre transform. Most textbooks on statistical mechanics do not tell its readers, however, the clear mathematical logic of all these formulae. But actually, Boltzmann’s 1877 paper [¹⁶], by counting the molecules with different kinetic energy in an ideal gas, had proceeded exactly the steps we took and derived the celebrated Boltzmann’s law, in the form in Eq. (11a).

Variational Bayesian method

The F(\beta ) obtained in Eq. (13) has a very important property: For any, arbitrary, normalized distribution \{z_k\},

(20){\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{z}_k\ln \left({\displaystyle \frac{z_k}{p_k{e}^{-\beta g_k}}}\right)\ge -\ln \left({\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^M}{p}_k{e}^{-\beta g_k}\right)=\beta F(\beta ).

In the variational Bayesian method for inference [¹⁷], one often knows a target, posterior distribution p_k{e}^{-\beta g_k} but computing its normalization factor is expensive. Eq. (20) shows that to obtain the target distribution, one can simply minimizes the left-hand-side of Eq. (20) among a set of possible \{z_k\}. This same idea had also been used by Gibbs in his variational method [¹⁸]: The free energy F(\beta ) of an equilibrium state is the minumum among all others through a virtual change of state.

Maximum entropy principle

The constrained optimization in Eq. (8) leading to distribu- tion in Eq. (11a) has also become the foundation of maximum entropy principle (MEP) championed by Jaynes [¹⁹], which has played a productive role in data science. The axiomatic nature of MEP [²⁰] and the role of conditional probability [²¹] have been elucidated.

The fundamental premises behind the large deviations principle (LDP) and the MEP are very different: Entropy, as a large deviation rate function, is used in the former to find the rare event that is the most probable, which is the only possible event in the limit: For an arbitrary set of n real values \{{\phi }_i\},

(21)(e^{-N_{{\phi }_1}}+e^{-N_{{\phi }_2}}+\cdots +e^{-N_{{\phi }_n}})\sim e^{-N_{{\phi }^{\ast }}}\mathrm{as}\hspace{0.17em}N\to \infty ,

where {\phi }^{\ast }=\min \{{\phi }_1,...,{\phi }_n\}. This is the same idea in choosing only the term with the largest eigenvalue among the terms in a linear eigenvalue decomposition, in the limit of infinite time or system’s size. In MEP, however, entropy function is used as a measure for “unbias”. Actually, according to LDP, the {\mathrm{x}}^{\ast } in (11a) is not a probability distribution, it is the most probable frequency among N i.i.d. samples. In MEP, it is interpreted as the least biased probability distribution with maximum uncertainty.

References

[1] Pence, C. H. (2011) “Describing our whole experience”: the statistical philosophies of W. F. R. Weldon and Karl Pearson. Stud. Hist. Philos. Biol. Biomed. Sci., 42, 475–485

[2] Chibbaro, S., Rondoni, L. and Vulpiani, A. (2014) Reductionism, Emergence and Levels of Reality. New York: Springer

[3] Anderson, P. W. (1972) More is different. Science, 177, 393–396

[4] Qian, H., Ao, P., Tu, Y. and Wang, J. (2016) A framework towards understanding mesoscopic phenomena: Emergent unpredictability, symmetry breaking and dynamics across scales. Chem. Phys. Lett., 665, 153–161

[5] Ge, H. and Qian, H. (2016) Mesoscopic kinetic basis of macroscopic chemical thermodynamics: A mathematical theory. Phys. Rev. E, 94, 052150

[6] Dembo, A. and Zeitouni, O. (1998) Large Deviations Techniques and Applications, 2nd ed. New York: Springer

[7] Kurtz, T. G. (1972) The relationship between stochastic and deterministic models for chemical reactions. J. Chem. Phys., 57, 2976–2978

[8] Liggett, T. M. (1985) Interacting Particle Systems. New York: Springer-Verlag

[9] Derrida, B. (1998) An exactly soluble nonequilibrium system: The asymmetric simple exclusion process. Phys. Rep., 301, 65–83

[10] Gang, H. (1986) Lyapunov function and stationary probability distribution. Zeit. Physik B: Cond. Matt. 65, 103–106

[11] Chang, R. and Goldsby, K. A. (2012) Chemistry, 11th ed. New York: McGraw-Hill

[12] Qian, H. (2019) Stochastic population kinetics and its underlying mathematicothermodynamics. In: The Dynamics of Biological Systems, Bianchi, A., Hillen, T., Lewis, M., Yi, Y. eds., pp. 149–188. Springer: New York

[13] Murray, J. D. (2011) Mathematical Biology II: Spatial Models and Biomedical Applications, 3rd ed. New York: Springer

[14] von Bertalanffy, L. (1950) The theory of open systems in physics and biology. Science, 111, 23–29

[15] Huang, K. (1963) Statistical Mechanics. New York: John Wiley & Sons

[16] Sharp, K. and Matschinsky, F. (2015) Translation of Ludwig Boltzmann’s paper “On the relationship between the second fundamental theorem of the mechanical theory of heat and probability calculations regarding the conditions for thermal equilibrium”. Entropy (Basel), 17, 1971–2009

[17] Ghahramani, Z. (2001) An introduction to hidden Markov models and Bayesian networks. Int. J. Pattern Recognit. Artif. Intell., 15, 9–42

[18] Pauli, W. (1973) Pauli Lectures on Physics: Thermodynamics and the Kinetic Theory of Gas. Cambridge: The MIT Press

[19] Jaynes, E. T. (2003) Probability Theory: The Logic of Science. London: Cambridge University Press

[20] Shore, J. E. and Johnson, R. W. (1980) Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy. IEEE Trans. Inf. Theory, 26, 26–37

[21] van Campenhout, J. M. and Cover, T. M. (1981) Maximum entropy and conditional probability. IEEE Trans. Inf. Theory, 27, 483–489

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature