Close to Optimal Cell Sensing Ensures the Robustness of Tissue Differentiation Process: The Avian Photoreceptor Mosaic Case

Abstract

The way that progenitor cell fate decisions and the associated environmental sensing are regulated to ensure the robustness of the spatial and temporal order in which cells are generated towards a fully differentiating tissue still remains elusive. Here, we investigate how cells regulate their sensing intensity and radius to guarantee the required thermodynamic robustness of a differentiated tissue. In particular, we are interested in finding the conditions where dedifferentiation at cell level is possible (microscopic reversibility), but tissue maintains its spatial order and differentiation integrity (macroscopic irreversibility). In order to tackle this, we exploit the recently postulated Least microEnvironmental Uncertainty Principle (LEUP) to develop a theory of stochastic thermodynamics for cell differentiation. To assess the predictive and explanatory power of our theory, we challenge it against the avian photoreceptor mosaic data. By calibrating a single parameter, the LEUP can predict the cone color spatial distribution in the avian retina and, at the same time, suggest that such a spatial pattern is associated with quasi-optimal cell sensing. By means of the stochastic thermodynamics formalism, we find out that thermodynamic robustness of differentiated tissues depends on cell metabolism and cell sensing properties. In turn, we calculate the limits of the cell sensing radius that ensure the robustness of differentiated tissue spatial order. Finally, we further constrain our model predictions to the avian photoreceptor mosaic.

1. Introduction

Decision-making is a process to identify important choices and responses which depends on some basic criteria [1]. Cell decision-making is a process where cells select a new state, such as cell fates or phenotypes, in response to their microenvironmental milieu. In this regard, pluripotent cell differentiation can be viewed as cell decision-making of inheritable fates. Typically, cells irreversibly acquire new fates by following a hierarchical lineage, where pluripotent stem cells in a proper microenvironment (stem cell niche) differentiate into, for example, bone, muscle, epithelial, and further specialized cells. The cell differentiation process encompasses the dramatic change of geometry, shape, gene expression inside the cell, etc. [2,3]. It is yet to be fully understood how the information available to pluripotent progenitors, including its intrinsically determined state and extrinsic microenvironmental signals, is encoded and processed by progenitors to generate different differentiated cell types.

In 2006, Takahashi and Yamanaka [4] discovered that almost any differentiated cell can be sent back in time to a state of pluripotency by expressing appropriate transcription factors (the Nobel Prize in Medicine 2012). Cell reprogramming can be externally induced via the delivery of transcription factors, naturally or in vitro [5,6]. Such a reversion of a differentiated to a pluripotent state is the idea behind cancer stem cells (CSC). The CSC theory proposes that, among all cancerous cells, a few act as stem cells that reproduce themselves and sustain cancer, much like normal stem cells that renew and sustain organs and tissues of the body. The process of somatic reprogramming using Yamanaka factors, where many of them are oncogenes, offers a glimpse into how cancer stem cells may originate. In particular, neurological cancers such as primary glioblastomas [7] and retinoblastomas [8] are resulting from dedifferentiation of the glial and photoreceptor (retinal) cells, respectively. In particular, retinoblastoma tumor cells lose their photoreceptorness and become malignant plastic cells, i.e., CSC [8]. Xu et al. also have shown that the cell of origin for retinoblastoma is a committed cone precursor—an almost terminally differentiated photoreceptor—that has lost an Rb gene, and not a pluripotent progenitor [9,10]. Here, we use the example of photoreceptor mosaic of avian retina to shed light on the following questions: (Q1) how do cell intrinsic dynamics and microenvironmental factors coordinate during development to produce organized tissues such as photoreceptor mosaics [11]; and (Q2) why is the differentiated mosaic so stable or how probable is the reversal of retinal tissue back to a pluripotent one, e.g., retinoblastoma?

The theory behind cell differentiation has been formalized in a metaphorical way by C. H. Waddington [12,13,14], which allows for developing a dynamical systems framework for modeling single-cell fate decisions [15,16,17]. Waddington has depicted the developmental process as a series of cell decisions that can be represented as bifurcations towards a differentiated state/phenotype. Practically, cell states, also called microstates, can be viewed as a vector of molecular expressions that are experimentally measured via high-throughput omics data, FACS, immunohistochemistry markers, etc. [14,18,19]. We note that such cellular microstates are technically different from the classical statistical mechanics definition of microstates. In general, differentiated states can be viewed as the fixed point(s) of microstate attractors [20,21]. Typically, these states can be associated with an appropriate probability distribution peaked around the fixed point, which is the deepest point of the valley in the Waddington potential. The situation is much more intricate in the case of pluripotent, stem cell-like states, where Waddington has also depicted them as attractors around a fixed point. However, this has been recently challenged by Furusawa and Kaneko [22], where they have shown that stem cell-like attractors can be viewed as limit-cycle attractors and not as stable fixed points. Biological observations of dynamic variability of single cells within pluripotent cell populations distinguish between pluripotency as a molecular state and pluripotency as a function, indicating that a pluripotent state is not unique but rather appears to be compatible with a wide variety of interchangeable molecular microstates (patterns of gene or protein expressions) [23,24]. In this view, pluripotent cells independently explore a variety of molecular expression states, and this phenotypic/state exploration transiently primes each individual cell to respond to a range of various differentiation-inducing stimuli, depending upon its instantaneous molecular state [7]. Finally, Waddington theory does not take into account cell sensing and the corresponding interactions that take place in a tissue. The aforementioned limitations of the Waddington approach require further theoretical development to answer the questions (Q1) and (Q2).

In order to answer (Q1), in this paper, we employ the recently proposed the Least microEnvironmental Uncertainty Principle (LEUP)—which is essentially a statistical mechanical theory for cell decision-making [25,26,27]—and apply it to the problem of cell differentiation. The LEUP is inspired by the theories of Bayesian brain hypothesis [28], the free-energy principle [29], and other dynamic Bayesian inference theories that try to explain human brain cognitive dynamics. Similar ideas have also been proposed in the influential work by Bialek [30]. Similar to these theories, the LEUP is based on the premise that cell internal molecular networks adapt to the sensed microenvironmental data and subsequently determine the relevant decisions. In turn, cells are encoding sensed information in genetic, epigenetic, translational, or transcriptional levels, where different timescales are related to each of these encoding levels, depending on the persistence of the microenvironmental stimuli. The reason to propose a theory such as the LEUP is the fact that the complexity of molecular networks does not allow us to know the exact involved dynamics, thus we use the LEUP as a kind of dynamic Bayesian inference to circumvent this complexity and to make predictions. The LEUP theory also implies a decrease in the local microenvironmental entropy of the cell decision-maker, which biologically translates into the actions of differentiating cells which lead to more organized tissues during development.

The central cellular process related to the LEUP is cell sensing. Cells can acquire knowledge about their microenvironment by various sensing mechanisms such as the binding of their receptors to diffusible ligands [31], pseudopodia extension [32], mechanosensing [33], proton-pump channels [34], gap junctions, etc. Cells can sense rapid changes of their milieus, where they mainly exploit two ways to decrease sensing errors: (1) by increasing the number of receptors or the responses of the downstream signaling pathways [35], and/or (2) by increasing sensing area [36]. The latter can span from the resting cell size to extensions via pseudopodia, blebbing, and other cell size regulation mechanisms. Pseudopodia or blebs can act as sensors, via a pressure sensing mechanism mediated by Piezo channels, allowing cells to decide when and where to migrate [37]. In this regard, we further specify our second question (Q2) as: can one calculate the limits of the cell sensing radius that ensure the robustness of differentiated tissue spatial order?

To answer (Q2), we develop a thermodynamic-like theory using the tools of stochastic thermodynamics for a generic cell differentiation process. The main biological assumption is that differentiated cells can reverse to an undifferentiated state, for instance, by the process of carcinogenesis, see also Refs. [8,9,38]; a recent review on the latest advances in research on the process of dedifferentiation both at cell and tissue levels can be found in [39]. Stochastic thermodynamics allows us to understand the conditions that, even though single-cell dedifferentiation is possible (microscopic reversibility), the system (tissue) is still able to be robust, and it maintains its spatial order and differentiation integrity (macroscopic irreversibility). Stochastic thermodynamics is a suitable tool for systems where small scale dynamics matters (e.g., soft matters, active matters, and biological systems); in such cases, the higher-order moments dominate [40,41,42,43,44]; for different experimental applications of stochastic thermodynamics, see Refs. [45,46]. In addition, formulating the laws of thermodynamics in the mesoscopic scale (specifically, at the level of trajectory) has also been investigated recently [47]. To describe the cell-level dedifferentiation process in the context of stochastic thermodynamics, we apply the fluctuation theorem, specifically the Crooks’ relation [48,49]. The fluctuation theorem, which can be considered as the heart of stochastic thermodynamics, initially had been derived to explain how irreversibility at the macroscopic level emerges from the underlying reversible dynamics and to estimate the probability of the violation of the second law of thermodynamics within a short amount of time for systems at small scales [50,51,52].

By combining elements of stochastic thermodynamics with the LEUP, we aspire to approach questions (Q1) and (Q2) within the context of the avian photoreceptor mosaic. This paper is organized as follows: in Section 2, we review the basic concepts of the LEUP and its connection to statistical mechanics. In Section 3, we demonstrate how optimal microenvironmental sensing is associated with differentiated tissue spatial configuration; therein, we examine our theory using the data obtained from the (avian) photoreceptor mosaic. We apply the fluctuation theorem to cell differentiation in Section 4 and demonstrate the thermodynamic robustness of this process in the case of the avian retina development. In Section 5, we show how cell sensing radius and total entropy production are related and determine the limits of the sensing radius in order for the tissue development to be robust. Finally, we conclude and discuss our results in Section 6.

2. Least microEnvironmental Uncertainty Principle

A cell collects information from its microenvironment and, based on that, it takes action (phenotypic decisions). In other words, a cell reacts to the environmental information by changing its own state. By adopting the notations of Ref. [25], we denote internal variables (such as gene expression, RNA molecules, metabolites, etc.) of the n-th cell as $x_n$ and external variables (such as chemical signals, ligands, homotypic or heterotypic interactions, cellular densities, translational proteins, etc.) of the n-th cell as $y_n$. The latter includes all the different extrinsic variables within the interaction radius of the n-th cell, i.e., $y_n=\left\{y\left(r\right):r\in (r_n,r_n+R]\right\}$, where $r_n$ is the position vector of the n-th cell, and R is the maximum microenvironmental sensing radius of the n-th cell. Cell sensing radius is a tunable variable controlled by the cell and is regulated by various biophysical mechanisms, as stated in the Introduction. This is precisely what R is modeling, where it has been considered as an intrinsic cell variable dictated by the LEUP dynamics. We also note that we have assumed that a sensing radius can be identified with the typical interaction radius considered in agent-based models. In Figure 1, we have shown the schematic diagram of the microenvironment of a differentiated cell as being composed of both pluripotent progenitors and differentiated cells.

We consider the behavior of cell as a Bayesian decision-maker that reacts to its microenvironment [25,53], such that

$$P\left(x_n\left(t+{\tau }_{\mathrm{decis}.}\right)\mid y_n\left(t\right)\right)=\frac{P\left(y_n\left(t\right)\mid x_n\left(t\right)\right)P\left(x_n\left(t\right)\right)}{P\left(y_n\left(t\right)\right)},$$

(1)

where ${\tau }_{\mathrm{decis}.}$ is the time needed for the cell to make a decision. The distribution $P\left(y_n\left(t\right)\mid x_n\left(t\right)\right)$ is the probability of microenvironmental information/data being collected by the cell at time t. In other words, it is the probability that the cell perceives all other cells, chemicals, and nutrients in its surroundings. The distribution $P\left(x_n\left(t\right)\right)$ is the prior probability of the cell current internal states. If these two distributions are multiplied and divided by the probability of external states $P\left(y_n\left(t\right)\right)$, then (1) implies that the resulting quantity is the posterior probability distribution of internal states $P\left(x_n\left(t+{\tau }_{\mathrm{decis}.}\right)\mid y_n\left(t\right)\right)$. The latter describes the most likely decision to be made by cells over the internal variables after processing the available information in the time period of ${\tau }_{\mathrm{decis}.}$.

As stated in the Introduction, we assume that the cell prior is continually being updated by the previous time step posterior, in the sense of Bayesian learning. We have also assumed a perfect transfer of the prior to the next time step posterior, although this process is very noisy. Now, by taking the logarithm of (1) and integrating over the joint probability distribution $P\left(x_n\left(t\right),y_n\left(t\right)\right)$, for small decision times, we obtain:

$$\frac{\partial S\left(x_n\mid y_n\right)}{\partial t}=\frac{1}{{\tau }_{\mathrm{decis}.}}\left[S\left(y_n\mid x_n\right)-S\left(y_n\right)\right]=-\frac{1}{{\tau }_{\mathrm{decis}.}}I\left(x_n,y_n\right),$$

(2)

where, for simplicity, we have dropped t as an argument and $I\left(x_n,y_n\right)$ is the mutual information. The above equation reaches equilibrium when the mutual information vanishes. Now, we make the first crucial assumption of our work (A1) that the cell decision time is much smaller than the asymmetric division time, i.e., ${\tau }_{\mathrm{decis}.}\ll {\tau }_{\mathrm{div}.}$. This can be justified as cell division, which is a prerequisite for differentiation that takes around ∼24 h; on the other hand, the relaxation timescale of the molecular networks responsible for deciding a new state is much shorter (∼1 h). This implies that the dynamics of internal variables $x_n$’s are much faster than the microenvironmental dynamics of $y_n$’s. According to this timescale separation, we can assume that variations of $S\left(y_n\mid x_n\right)$ and $S\left(y_n\right)$ belong to the slow manifold of the system. In turn, we observe that in order for Equation (2) to reach an equilibrium, the microenvironmental entropy sensed by cells $S\left(y_n\mid x_n\right)$ should be inevitably a decreasing quantity with respect to time. This fact can be interpreted as the entropy of cell sensing distribution needing to become more focused and more independent of the microenvironment as time goes on, since this information—microenvironmental data—has already been encoded in the cell prior. The latter is valid within biophysical context because of the following reasons: first, the collection of microenvironmental data, that is, the precise evaluation of $P\left(y_n\mid x_n\right)$ is energetically expensive, thus it is not favorable as cells would like to be energetically efficient; and, second, fully differentiated cells are perfectly adapted to their microenvironment and the corresponding fluctuations imply an optimal prior distribution $P\left(x_n\right)$.

Using the timescale separation of (A1), we formulate a variational problem concerning the maximization of entropy of cell microstates $S\left(x_n\right)$ in order to find an optimal prior $P\left(x_n\right)$, which is constrained by the current local microenvironmental entropy $\overline{S}\left(y_n\mid x_n\right)$ that corresponds to a random sample of size R of the fully differentiated tissue of interest. We note that, due to (A1), the cell “solves” the entropy maximization by sensing an almost time-invariant microenvironmental entropy. We consider the definition of entropy in terms of probability distributions as $S=-{\Sigma }_i{P}_{i}ln{P}_i$ (Shannon entropy), where, in the case of a continuous set of states, the summation symbol is to be replaced by an integral. Thus, our variational problem reads

$$\begin{array}{cc}\hfill \frac{\delta }{\delta P\left(x_n\right)}& \{S\left(x_n\right)-\beta \left[\int P\left(x_n\right)S\left(y_n\mid x_n\right)\mathrm{d}{x}_n-\overline{S}\left(y_n\mid x_n\right)\right]\hfill \\ \hfill & -\lambda \left[\int P\left(x_n\right)\mathrm{d}{x}_n-1\right]\}=0,\hfill \end{array}$$

(3)

where $\delta /\delta P\left(x_n\right)$ is the functional derivative. We have introduced two Lagrange multipliers, namely, $\beta$ and $\lambda$. The parameter $\beta$ is related to the current local microenvironmental entropy $\overline{S}\left(y_n\mid x_n\right)$ and $\lambda$ is a constraint which preserves the normalization of probability. We note that other knowledge about the system, such as an explicit model in terms of statistical observables, can be incorporated as additional constraints. The solution of Equation (3) then reads

$$P_1\left(x_n\right)=\frac{e^{-\beta S\left(y_n\mid x_n\right)}}{Z_1},$$

(4)

where $Z_1\left(\beta \right)=\int e^{-\beta S\left(y_n\mid x_n\right)}\mathrm{d}{x}_n$ is a normalization constant. We note that the parameter $\beta$ quantifies the compliance of the cell to the LEUP, which is related to the ability of cells to sense their microenvironment and convert the perceived effects into biological signals. In the next section, we elaborate on the relation of $\beta$ and microenvironmental sensing.

Now, we want to generalize (4) to a biologically relevant scenario of cell differentiation. It is well known that differentiation takes place during asymmetric divisions of pluripotent progenitor cells [54]. Let us assume that $\mu \propto {\tau }_{\mathrm{div}.}$ is the asymmetric division probability of a pluripotent cell. In a mesoscopic microenvironmental ensemble of $N\left(y_n\mid x_n\right)$ cells around the n-th cell, we want to calculate the probability of the central n-th cell to change its phenotype. This probability reads as

$$\begin{array}{cc}\hfill P\left(x_n\right)=& \mathrm{Prob}.\left\{\mathrm{to}\mathrm{select}\mathrm{the}n-\mathrm{th}\mathrm{cell}\mathrm{out}\mathrm{of}N\mathrm{cells}\right\}\hfill \\ \hfill & \times \mathrm{Prob}.\left\{\mathrm{the}n-\mathrm{th}\mathrm{cell}\mathrm{divides}\mathrm{asymmetrically}\mathrm{once}\right\}\hfill \\ \hfill & \times \mathrm{Prob}.\left\{\mathrm{the}n-\mathrm{th}\mathrm{cell}\mathrm{decides}\mathrm{over}\mathrm{its}\mathrm{phenotype}\right\}.\hfill \end{array}$$

(5)

The latter probability happens according to the LEUP, and it is the same as (4). Now, we assume that the probability of asymmetric divisions within an ensemble of N cells follows a Poisson distribution Pois$\left[\mu N\left(y_n\mid x_n\right)\right]$, which is a limit of the Binomial distribution B$[N\left(y_n\mid x_n\right),\mu ]$ for a small proliferation probability $\mu$. Putting all these together, we can write:

$$P\left(x_n\right)\propto \frac{1}{N\left(y_n\mid x_n\right)}\times N\left(y_n\mid x_n\right)e^{-\mu N\left(y_n\mid x_n\right)}\times e^{-\beta S\left(y_n\mid x_n\right)}=\frac{e^{-\beta S\left(y_n\mid x_n\right)-\mu N\left(y_n\mid x_n\right)}}{Z\left(\beta ,\mu \right)},$$

(6)

where the normalization factor of $P\left(x_n\right)$ is defined as $Z\left(\beta ,\mu \right)=\int e^{-\beta S\left(y_n\mid x_n\right)-\mu N\left(y_n\mid x_n\right)}\mathrm{d}{x}_n$.

We illustrate the above with a concrete example: if the microenvironmental probability distribution sensed by the n-th cell follows a Gaussian distribution with the variance ${\sigma }_n^2$, assuming that $y_n$ is scalar, and hence $S\left(y_n\mid x_n\right)=(1/2)ln\left[2\pi e{\sigma }_n^2\left(y_n\mid x_n\right)\right]$, and, by considering $\mu =0$, then (6) reduces to

$$P\left(x_n\right)=\frac{{\sigma }_n^{-\beta }\left(y_n\mid x_n\right)}{{\sum }_j{\sigma }_j^{-\beta }\left(y_j\mid x_j\right)}.$$

(7)

Now, we are in a position to establish an explicit formula for cell internal entropy. To this end, by exploiting (6), we obtain

$$S\left(x_n\right)=-\int P\left(x_n\right)lnP\left(x_n\right)\mathrm{d}{x}_n=\beta {\langle S\left(y_n\mid x_n\right)\rangle }_{x_n}+\mu {\langle N\left(y_n\mid x_n\right)\rangle }_{x_n}+lnZ,$$

(8)

where $\langle \dots \rangle$ denotes the expectation value. By making usage of the notions which are developed within the context of thermodynamics, we can also define phenotypic internal energy sensed by the cell as

$$U_{\mathrm{int}.}\left(x_n,y_n\right)={\langle S\left(y_n\mid x_n\right)\rangle }_{x_n}+\frac{\mu }{\beta }{\langle N\left(y_n\mid x_n\right)\rangle }_{x_n}.$$

(9)

We note that, for $\mu =0$, the internal energy is the same as entropy for a particular realization. Thus, (6) resembles the Boltzmann distribution.

Interestingly, since the microenvironmental entropy $S\left(y_n\mid x_n\right)$ acts as an effective internal energy, it is expected to be a decreasing quantity in time. This result is in agreement with the discussion which follows (2). Thus, the LEUP is consistent with its premises.

The LEUP and Known Statistical Results

The steady-state distribution of (7) which is derived by the LEUP, under the assumption of Gaussian distributed microenvironment and $\mu =0$, can be associated with different established statistical results. Here, we distinguish several cases for specific values of $\beta$, which have particular meanings in statistics. For instance, when $\beta =1$, we reproduce distributions proportional to the well-known Jeffreys prior, see also Ref. [25]. Jeffreys prior is known as the most typical uninformative prior used in Bayesian inference [55].

Interestingly, if we assume that $\beta =2$, then we recover the so-called minimum variance estimator when fusing multiple scalar estimates. Following [56], let us assume extrinsic variables $y_1,y_2,\dots ,y_n$, such as ligand concentrations, cell densities, etc. that are pairwise statistically uncorrelated and are normally distributed, i.e., $y_i\sim \mathcal{N}({\mu }_i,{\sigma }_i^2\left(x_i\right))$, where $1\le i\le n$. In turn, we assume intrinsic variables $x_1,x_2,\dots ,x_n$ that correspond to cellular sensors and a downstream process of the aforementioned microenvironmental variables. We define the average sensed microenvironment as the average of the extrinsic signals $y_i$’s weighted by the distribution of the cell sensors $x_i$’s, i.e., $Z(y_1,y_2,\dots ,y_n)={\sum }_{j=1}^n{P}_j{y}_j$, where ${\sum }_{j=1}^n{P}_j=1$. The distribution of internal variables $P_i$ that minimizes the variance, aka noise of the sensed microenvironment, is given by the following formula:

$$P_i=\frac{{\sigma }_i^{-2}\left(x_i\right)}{{\sum }_j{\sigma }_j^{-2}\left(x_j\right)},$$

(10)

which has the same form as the steady state of the LEUP (7), for $\beta =2$. In the following, we connect the latter result with a specific cell sensing scenario.

3. Minimal Cell Sensing Noise Is Associated with Differentiated Tissue Spatial Configuration

In this section, first we show how the LEUP is connected to cell sensing with respect to receptor–ligand binding. Then, we apply our theory to avian photoreceptor mosaic and fit the parameter $\beta$ to recover the photoreceptor percentages in retina.

3.1. How Does the LEUP Relate to Cell Sensing?

In this subsection, we show the relevance of the LEUP within the context of cell sensing mechanisms. Here, we focus on the receptor–ligand sensing apparatus of cell, where complexes are formed at the cell membrane and subsequently internalized via endocytosis [31]. This sort of sensing mechanism is also relevant in the case of photoreceptors [57]. Receptors act as sensors for specific microenvironmental molecules (ligands) that bind together with a certain affinity to form complexes. Often, different ligands may bind to the same receptor, for instance, the Notch-Delta-Jagged system, where the Notch receptor can bind to either Delta or Jagged molecules [58,59]. The sensed information can be quantified by the concentration of internalized complex molecules. For simplicity, we denote x as receptor concentration at the cell membrane and $y_1$ and $y_2$ as the corresponding ligands concentrations, which we assume are statistically independent. Moreover, the ligands’ concentrations are consumed by cells, since they bind to receptors and therefore depend on the concentration x (we omit this dependence for the sake of simplicity). Typical dynamics of such systems read as

$$\frac{dc}{dt}=k_{1}x{y}_1+k_{2}x{y}_2-dc.$$

(11)

Under the assumption of fast decay rates $d\gg 1$, which is consistent with (A1), the system behaves as in a steady state, that is,

$$c_{\mathrm{eq}.}=\frac{k_1}{d}x{y}_1+\frac{k_2}{d}x{y}_2.$$

(12)

The terms, $(k_i/d)x,i=1,2$, define the percentage of receptors bound to ligands $y_i$’s. In the context of the LEUP, $(k_i/d)x$ corresponds to $P(x\mid y_i)$, which, under the assumption of Gaussian microenvironment—in a steady state—is:

$$P_i=\frac{{\sigma }_i^{-\beta }\left(x_i\right)}{{\sigma }_1^{-\beta }\left(x_1\right)+{\sigma }_2^{-\beta }\left(x_2\right)},i=1,2,$$

(13)

which for $\beta =2$ coincides with (10). Taking all the above together, we translate the receptor–ligand cell sensing system as a linear combination of complex formation estimates $y_i$’s, where the LEUP probabilities $P_i$’s are the corresponding coefficients/proportions of complexes bound to ligands $y_i$’s, that is,

$$c_{\mathrm{eq}.}=\frac{{\sigma }_1^{-\beta }\left(x_1\right)}{{\sigma }_1^{-\beta }\left(x_1\right)+{\sigma }_2^{-\beta }\left(x_2\right)}{y}_1+\frac{{\sigma }_2^{-\beta }\left(x_2\right)}{{\sigma }_1^{-\beta }\left(x_1\right)+{\sigma }_2^{-\beta }\left(x_2\right)}{y}_2,$$

(14)

which coincides with the definition of $Z(y_1,y_2,\dots ,y_n)$, for $n=2$, defined in the previous section.

At this point, we want to do a sanity test for our results. For simplicity, we consider a single type of ligand concentration y that binds to a receptor of concentration x. Then, the complex formation dynamics are described as

$$\frac{dc}{dt}=k_{+}xy-dc.$$

(15)

At the equilibrium, the steady-state complex concentration $c_{\mathrm{eq}.}$ reads

$$c_{\mathrm{eq}.}=\frac{k_{+}}{d}xy.$$

(16)

As before, we define the probability $P_x=(k_{+}/d)x$ as the proportion of binding receptors. Then, from (16), we observe that:

$$P_x\propto \frac{1}{y}.$$

(17)

The ligand concentration, $y=\langle Y\rangle$, is the expected value of the number of ligand molecules within the cell volume, which is denoted as Y. Using arguments similar to that of Berg and Purcell in their seminal paper [35], the ligand molecules are diffusing and therefore we can assume that Y follows a Poisson distribution, thus ${\sigma }_Y^2=y$, where the ligand variance is denoted as ${\sigma }_Y^2$. Now, by combining the latter with (17), we recover our non-normalized the LEUP result for $\beta =2$:

$$P_x\propto {\sigma }_y^{-2}\left(x\right).$$

(18)

This specific value of $\beta$ complies with the perfect monitoring/sensing assumptions (for details, see [35]). We note that, in (18), we have explicitly denoted the dependence of diffusible ligand concentration on the receptor concentration x. As a side remark, if we change the receptor–ligand binding term by introducing, for instance, finite number of receptors or covalent bonds, then the parameter $\beta$ will be modified.

The above result is pivotal since it connects the cell fate decision-making with optimal microenvironmental sensing in terms of minimization of sensing noise, which results in a specific spatial phenotypic distribution. In the following, we explore the validity of these results in the case of the avian photoreceptor mosaic.

3.2. Application: Predicting the Cone Color Distribution in the Avian Photoreceptor Mosaic

In this subsection, we focus on a particular multicellular system, namely, the photoreceptor cone mosaic of the avian retina. Retina plays a vital role in the visual perception of vertebrates allowing for light sampling and signal transfer to the corresponding optic nerves. The photoreceptor layer is responsible for the efficient light collection and color encoding of the surrounding structures. The avian retina possesses one of the most sophisticated cone photoreceptor systems among vertebrates [60]. Birds have five types of cones including four single cones, namely, green, red, blue, and violet, and a double cone, whereas, in humans, there are only three: red, green, and blue. In the case of birds, these cones are ordered on the retina in a non-perfect hexagonal manner. The five avian cone types exist as five independent spatial mosaics that are all embedded within a single monolayered epithelium along with a population of rod photoreceptors. In Ref. [60], the spatial organization of these photoreceptor mosaics has been analyzed, where the authors have found that the five cone types are present in the characteristic ratios at the latest stage of development, as is illustrated in Figure 2.

Our goal is to reproduce the ex vivo observed cone distributions of [60] by using the aforementioned LEUP-driven results. In particular, we intend to fit the parameter $\beta$ in formula (7) using the photoreceptors data. Therefore, let us consider that we have $n=5$ different cone states $x_i$’s that correspond to the expressions of different opsins, $i=g,r,b,v,\delta$ (g stands for green, r for red, b for blue, v for violet, and $\delta$ for double). Moreover, we define the corresponding cone occurrence probabilities $P_i=P\left(x_i\right)$ as cone photoreceptors’ percentages in retina. In turn, we denote ${\sigma }_i$ as the standard deviation for a given spatially local neighborhood related to the Nearest Neighbor Distribution (NND). In Ref. [60], the authors have exploited the Delaunay triangulation in order to calculate the NND for each cone and they have reported the first and second moments of it. In order to compare the LEUP color distributions against the ex vivo ones, we employ the Kullback–Leibler divergence. In turn, we numerically search for the $\beta$’s which appropriately fit the physiologically observed photoreceptor fate ratios. In the left panel of Figure 3, we have depicted the Kullback–Leibler divergence, i.e., $D_{KL}={\sum }_i{P}_{i}ln(P_i/Q_i)$, as a function of $\beta$, where i denotes different cone cells, and $P_i$ and $Q_i$ correspond to the probabilities obtained from the experimental data and the LEUP, respectively. The plot demonstrates that $D_{KL}$ has a minimum of $\approx 0.004$ at $\beta \approx 1.754$. We have compared the ex vivo observed cone distributions to the LEUP, in the right panel of Figure 3. Interestingly, the $\beta$ values that provide the best fit are in the range $1.378\le \beta \le 2.136$, where the values of the Kullback–Leibler divergence are on the order of magnitude ${10}^{-3}$. Thus, it demonstrates that the values of $\beta \approx 2$ explain the observed data, implying that the avian cone mosaic distribution is associated with an optimal cell sensing process that minimizes the corresponding noise.

4. Thermodynamic Robustness of Differentiated Tissues: The Fluctuation Theorem

In this section, we determine the thermodynamic constraints of two coarse-grained cell states that correspond to pluripotent (s) and differentiated (d) states. Then, we apply the obtained results to the particular case of the avian cone cells’ differentiation and the formation of photoreceptor mosaics. To this end, first we show how microstates (internal variables) are related to the microenvironmental information and heat transfer (cell metabolism).

In this context, a cellular microstate corresponds to a cell phenotype that lives in a tissue, which could be gene expression, RNA molecules, receptor distribution, etc. In other words, microstate gives information about the internal states of cell. We label these internal variables as $x_s$ and $x_d$ corresponding to pluripotent and differentiated cells, respectively. We define a cellular macrostate as a statistical observable (e.g., average) of a cell microenvironment that involves multiple cells of different phenotypes. Macrostates contain information about external variables which are labeled as $y_s$ and $y_d$. The macrostate $y_s$ is assumed to describe a microenvironment of pluripotent progenitor or stem cells characterized by the microstate $x_s$; the macrostate $y_d$ accordingly describes a microenvironment of differentiated cells characterized by the microstate $x_d$. We denote the number of pluripotent cells neighboring a cell with microstate $x_s$ as $N(y_s\mid x_s)$ and the number of differentiated cells neighboring a cell with microstate $x_d$ as $N(y_d\mid x_d)$. The total number of pluripotent and differentiated cells inside system are denoted as $N\left(s\right)$ and $N\left(d\right)$, respectively.

Now, based on (6), we can write the probability of cell to be in the microstate $x_s$ with the corresponding macrostate s as

$$P\left(x_s\right)=\frac{e^{-{\beta }_{s}S\left(y_s\mid x_s\right)-{\mu }_{s}N\left(y_s\mid x_s\right)}}{Z_1},$$

(19)

and the probability of being in the microstate $x_d$ with the corresponding macrostate d as

$$P\left(x_d\right)=\frac{e^{-{\beta }_{d}S\left(y_d\mid x_d\right)-{\mu }_{d}N\left(y_d\mid x_d\right)}}{Z_2}.$$

(20)

It is known that, for a system which is coupled to a set of heat baths and is in a time-symmetrically driven non-equilibrium state, the Crooks’ theorem is applicable [48,49]. In this case, the dynamics follow Brownian motion, i.e., there are no extreme jumps in the system state. For any trajectory to be initially at $x_s\left(0\right)$ and going through microstates $x_s\left(t\right)$ over time $\tau$, the Crooks’ theorem implies that [49],

$${\beta }^{{}^{\prime }}\Delta Q=ln\left[\frac{w\left[x_s\left(t\right)\right]}{w\left[x_s(\tau -t)\right]}\right],$$

(21)

where ${\beta }^{{}^{\prime }}\equiv 1/T$, for which T is the temperature of heat bath, $\Delta Q$ is the total heat released into bath over the course of $x_s\left(t\right)$, and $w\left[x_s\left(t\right)\right]$ is the probability of trajectory $x_s\left(t\right)$. Equation (21) demonstrates that, when there is a forward state change, the system loses heat to reservoir and, in the case of a time-reversed path, there is a heat gain from reservoir; this, in turn, implies that a forward trajectory is more probable than a time-reversed one; thus, Equation (21) substantiates a relation between heat and irreversibility at a microscopic (cell) level. This conclusion is in agreement with the ideas presented in [7]. The authors of the paper mentioned have considered differentiation as a series of reversible transitions through many microstates (we shall use this notion when discussing differentiation at a tissue scale in the following), where stem cells exhibit reversible oscillations until an attractor drives them towards a differentiated state. Within this picture, dedifferentiation is more likely to occur only on a small scale with a low probability, as a series of microstate transitions. We also note that, as is demonstrated in Ref. [49], Equation (21) is also valid for other kinds of steady-state probability distributions besides the classical Boltzmann distribution. In our case, such Brownian jumps are interpreted as changes in phenotypes and the associated heat losses are assumed to be due to potential changes in cell metabolism. Another important heat loss contribution comes from the cell divisions, which are required in the differentiation process [62]. However, there are other minor heat loss sources that are disregarded since they act on shorter time scales, such as physical friction, changes in the cytoskeleton, etc.

Our goal here is to establish a microreversibility relation for a general differentiation process (the details of derivation can be found in Appendix A). To this end, first, we fix the starting point of a differentiation trajectory as $x_s$, which is a microstate realization that belongs to a pluripotent cell attractor/set, and the ending point as $x_d$ that is a member of the differentiated phenotypes set (see Figure A1 in Appendix A). The Crooks’ theorem in (21) is valid for a single differentiation trajectory that connects a single pluripotent state to a particular differentiation state realization out of many. Therefore, first, we average over all possible paths that connect a pluripotent and a differentiation state, and then we average over combinations of dedifferentiation paths that connect any pluripotent to any differentiation states, resulting in:

$$\frac{w(x_s\to x_d)}{w(x_d\to x_s)}={\langle exp[{\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }]\rangle }_{x_d\to x_s},$$

(22)

where $w(x_s\to x_d)$ is the transition probability for which the system is found to be in the microstate $x_d$ at time $\tau$, given that the system was initially in the microstate $x_s$. The averaged version of the Crooks’ theorem in (22) has an important implication since it demonstrates that paths leading to differentiation and also dedifferentiation are many more than one.

Now that we have a microscopic relation of (22) at our disposal, we can study macroscopic (tissue) consequences of that. At this point, we make the second vital assumption of our study (A2) that the tissue dynamics follow a Markov process. The probabilistic description of macroscopic states of two distinct cell types, from which one can understand the phenomenon of irreversibility at a macroscopic level, can be constructed as [49]:

$$W(s\to d)={\int }_d\mathrm{d}{x}_d{\int }_s\mathrm{d}{x}_{s}P(x_s\mid s)w(x_s\to x_d),$$

(23)

and

$$W(d\to s)={\int }_s\mathrm{d}{x}_s{\int }_d\mathrm{d}{x}_{d}P(x_d\mid d)w(x_d\to x_s),$$

(24)

where $P(x_s\mid s)$ is the probability that the system is in the microstate $x_s$, given that it is observed in the macrostate s, and $w(x_s\to x_d)$ is defined as before, where we have omitted $\tau$ for the notational convenience. The transition probability $W(s\to d)$ in (23) implies the likelihood that the cell can be observed in the macrostate d while it was initially prepared in the macrostate s; accordingly, (24) is understood in the same fashion, i.e., the likelihood that a microenvironment is prepared in the state d to satisfy the microenvironment s after another time interval $\tau$. These processes are illustrated in Figure 4.

By taking the ratio of (23) and (24), we obtain

$$\frac{W(d\to s)}{W(s\to d)}=\frac{{\int }_s\mathrm{d}{x}_s{\int }_d\mathrm{d}{x}_d[P(x_d\mid d)/P\left(x_d\right)]P\left(x_d\right)w(x_d\to x_s)}{{\int }_d\mathrm{d}{x}_d{\int }_s\mathrm{d}{x}_s[P(x_s\mid s)/P\left(x_s\right)]P\left(x_s\right)w(x_s\to x_d)},$$

(25)

where we have multiplied and divided the numerator and the denominator by $P\left(x_d\right)$ and $P\left(x_s\right)$, respectively. We define the pointwise mutual information for individual trajectories as $i_1=ln[P(x_s\mid s)/P\left(x_s\right)]$ and $i_2=ln[P(x_d\mid d)/P\left(x_d\right)]$, and then, by taking these definitions into account and replacing $P\left(x_d\right)$ by its corresponding relation in (20) and (25) reduces to

$$\begin{array}{cc}\hfill \frac{W(d\to s)}{W(s\to d)}& =\frac{{\int }_s\mathrm{d}{x}_s{\int }_d\mathrm{d}{x}_d{e}^{i_2}[e^{-{\beta }_{d}S\left(y_d\mid x_d\right)-{\mu }_{d}N\left(y_d\mid x_d\right)}/Z_2]w(x_d\to x_s)}{{\int }_d\mathrm{d}{x}_d{\int }_s\mathrm{d}{x}_s{e}^{i_1}P\left(x_s\right)w(x_s\to x_d)}\hfill \\ \hfill & =\frac{{\int }_s\mathrm{d}{x}_s{\int }_d\mathrm{d}{x}_d{e}^G{e}^{i_2+i_1-i_1}(Z_1/Z_2)[e^{-{\beta }_{s}S\left(y_s\mid x_s\right)-{\mu }_{s}N\left(y_s\mid x_s\right)}/Z_1]w(x_d\to x_s)}{{\int }_d\mathrm{d}{x}_d{\int }_s\mathrm{d}{x}_s{e}^{i_1}P\left(x_s\right)w(x_s\to x_d)},\hfill \end{array}$$

(26)

where

$$G\equiv -[{\beta }_{d}S\left(y_d\mid x_d\right)-{\beta }_{s}S\left(y_s\mid x_s\right)]-[{\mu }_{d}N\left(y_d\mid x_d\right)-{\mu }_{s}N\left(y_s\mid x_s\right)].$$

By exploiting (19) and (22), we can rewrite (26) as

$$\frac{W(d\to s)}{W(s\to d)}=\frac{{\int }_s\mathrm{d}{x}_s{\int }_d\mathrm{d}{x}_d{e}^G{e}^{i_1}{e}^{\Delta i}(Z_1/Z_2)P\left(x_s\right){\langle e^{-{\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }}\rangle }_{x_d\to x_s}w(x_s\to x_d)}{{\int }_d\mathrm{d}{x}_d{\int }_s\mathrm{d}{x}_s{e}^{i_1}P\left(x_s\right)w(x_s\to x_d)},$$

(27)

where $\Delta i\equiv i_2-i_1$.

Now, (27) can be expressed in terms of the average over all trajectories from the ensemble of microstates $x_s$ which correspond to the macrostate s to the ensemble of microstates $x_d$ which correspond to the macrostate d while each path is weighted by its probability as

$$\begin{array}{cc}\hfill \frac{W(d\to s)}{W(s\to d)}=\langle & {\langle exp[-{\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }]\rangle }_{x_d\to x_s}exp[-\Delta S_{\mathrm{LEUP}}]\hfill \\ \hfill & \times exp[-\Delta N]exp[\Delta i]exp[ln(Z_1/Z_2)]{\rangle }_{s\to d},\hfill \end{array}$$

(28)

where

$$\begin{array}{cc}\hfill & \Delta S_{\mathrm{LEUP}}\equiv {\beta }_{d}S(y_d\mid x_d)-{\beta }_{s}S(y_s\mid x_s),\hfill \\ \hfill & \Delta N\equiv {\mu }_{d}N(y_d\mid x_d)-{\mu }_{s}N(y_s\mid x_s).\hfill \end{array}$$

(29)

By rearranging (28) as

$$\begin{array}{cc}\hfill \langle & {\langle exp[-{\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }]\rangle }_{x_d\to x_s}exp[-\Delta S_{\mathrm{LEUP}}]exp[-\Delta N]exp[\Delta i]\hfill \\ \hfill & \times exp[ln(Z_1/Z_2)]exp[ln\{W(s\to d)/W(d\to s)\}]{\rangle }_{s\to d}=1,\hfill \end{array}$$

(30)

and using the Jensen’s inequality, i.e., $exp[\langle X\rangle ]\le \langle exp\left[X\right]\rangle$, and the fact that $e^x\ge 1+x$, we arrive at

$$\begin{array}{cc}\hfill {\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}& +{\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}+{\langle \Delta N\rangle }_{s\to d}-{\langle \Delta i\rangle }_{s\to d}\hfill \\ \hfill & -ln\left[\frac{Z_1}{Z_2}\right]+ln\left[\frac{W(d\to s)}{W(s\to d)}\right]\ge 0.\hfill \end{array}$$

(31)

We note that, if s and d correspond to the same identical classes, i.e., $W(d\to s)=W(s\to d)$, then the last term of (31) vanishes. Now, by defining the remaining terms as total entropy production, that is,

$$f\equiv {\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}+{\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}+{\langle \Delta N\rangle }_{s\to d}-{\langle \Delta i\rangle }_{s\to d}-ln\left[\frac{Z_1}{Z_2}\right],$$

(32)

then (31) implies that total entropy production is always non-negative. In other words, (31) can be considered as a generalization of the second law of thermodynamics (see also Ref. [49]).

The above inequality is reformulation of the fluctuation theorem for tissue differentiation. It relates thermodynamic properties of the system to the LEUP; and, as is obtained under general assumptions, it is applicable to a general cell and tissue differentiation process. In the next subsection, we show how (31) implies the thermodynamic robustness of cell differentiation, in the particular case of the avian cone photoreceptor mosaics formation.

Application: Differentiated Photoreceptor Mosaics Are Thermodynamically Robust

In this subsection, we use the fluctuation theorem, culminated in (31), to demonstrate the robustness of the avian photoreceptor mosaics. By defining $s\to d$ as a pluripotent tissue that differentiates to a more specialized one and denoting its corresponding forward transition probability as $p_f$, i.e., $W(s\to d)=p_f$, and $d\to s$ as differentiated tissue dedifferentiates to the pluripotent one with its backward transition probability of $p_b$, i.e., $W(d\to s)=p_b$, we can rewrite (31) in terms of $p_f$ and $p_b$ as

$$\frac{p_f}{p_b}\le exp\left[{\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}+{\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}+{\langle \Delta N\rangle }_{s\to d}-{\langle \Delta i\rangle }_{s\to d}-ln\left[\frac{Z_1}{Z_2}\right]\right],$$

(33)

where we have already defined the exponent of the exponential on the right-hand side of (33) as total entropy production, see (32).

In order to obtain a thermodynamic constraint which ensures the robustness of differentiated tissue, first we assume that there exists a maximum forward transition probability from s to d in such a way that

$$\frac{{p_f}_{\max }}{p_b}=exp\left[{\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}+{\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}+{\langle \Delta N\rangle }_{s\to d}-{\langle \Delta i\rangle }_{s\to d}-ln\left[\frac{Z_1}{Z_2}\right]\right].$$

(34)

To simplify (34) more, we note that, based on (8), we can write

$$\begin{array}{cc}\hfill & S\left(x_s\right)={\beta }_s{\langle S\left(y_s\mid x_s\right)\rangle }_{x_s}+{\mu }_s{\langle N\left(y_s\mid x_s\right)\rangle }_{x_s}+ln{Z}_1,\hfill \\ \hfill & S\left(x_d\right)={\beta }_d{\langle S\left(y_d\mid x_d\right)\rangle }_{x_d}+{\mu }_d{\langle N\left(y_d\mid x_d\right)\rangle }_{x_d}+ln{Z}_2,\hfill \end{array}$$

(35)

and then, by subtracting these two and taking the average over all trajectories from s to d, we obtain

$${\langle \Delta S\rangle }_{s\to d}\equiv {\langle S\left(x_d\right)-S\left(x_s\right)\rangle }_{s\to d}={\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}+{\langle \Delta N\rangle }_{s\to d}-ln\left[\frac{Z_1}{Z_2}\right],$$

(36)

where we have used the definitions of $\Delta S_{\mathrm{LEUP}}$ and $\Delta N$ as given in (29). Now, (34) can be written as

$$\frac{{p_f}_{\max }}{p_b}=exp\left[{\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}+{\langle \Delta S\rangle }_{s\to d}-{\langle \Delta i\rangle }_{s\to d}\right].$$

(37)

The immediate implication of (37) is that differentiated tissue is thermodynamically robust if

$${\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}+{\langle \Delta S\rangle }_{s\to d}>{\langle \Delta i\rangle }_{s\to d}.$$

(38)

In the following, we present the avian cone cells’ differentiation and the corresponding photoreceptor mosaics formation as an example of which the above inequality is satisfied. In other words, we demonstrate the robustness of the avian retina development and the irreversibility of the time arrow in this particular process.

First, we note that, as heat dissipation (the first term in the left-hand side of (38)) depends on metabolic pathways, this implies the crucial role of cell metabolism in the process of cell differentiation. In addition, consumption of glucose depends upon cell types. Progenitor or pluripotent stem cells use glucose as the primary metabolites for anaerobic glycolysis (fermentation) pathway: Glucose + 2 ADP + 2 Phosphate → 2 Lactate + 2 H${}^{+}$ + 2 ATP, where, as a result of the breakdown of glucose to lactic acid, the amount of energy around $109.4$ kJ/mol is released [63]. Differentiated cells use glucose to produce carbon dioxide and water by using aerobic (respiration) reaction: Glucose + 6O${}_2$ + 36ADP + 36Phosphate → 6CO${}_2$ + 6H${}_2$O + 36ATP, where the released energy is around 2820 kJ/mol [63]. Thus, ${\langle \Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}$, which is the total heat released during the journey from s to d is $2929.4$ kJ/mol. As these reactions have taken place in $T=310$ K, we have: ${\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}\approx 9.450$, where we have dropped its dimension as we have set the Boltzmann constant to be 1, in this paper.

Now, in order to investigate (38) for the avian cone cells’ differentiation, we need to calculate the value of ${\langle \Delta S\rangle }_{s\to d}$ for this particular process. Within the context of photoreceptor cells, differentiated states are given by an array of five opsin expressions as $x=(x_g,x_r,x_b,x_v,x_{\delta })$ and the corresponding probability distribution of $P\left(x\right)=(P_g,P_r,P_b,P_v,P_{\delta })$. At this point, we assume—(A3)—that pluripotent state corresponds to a state which is compatible to an equiprobable distribution of opsins, that is, $P\left(x_s\right)=(1/5,1/5,1/5,1/5,1/5)$. As is discussed in the Introduction, pluripotent states are not considered as attractive fixed point(s) but rather as oscillating attractors that explore the cellular state space [22]. Differentiated states $x_d$’s are all the states that their $P\left(x_d\right)$’s are equal to unimodal distributions centered at certain colors. Kram et al. [60] have reported different color percentages of the green, red, blue, violet, and double cone cells inside the avian retina (see Supporting Information of the mentioned reference). By assuming these numbers as the occurrence probabilities of the corresponding cones, we can write: $P_g\approx 0.204$, $P_r\approx 0.160$, $P_b\approx 0.133$, $P_v\approx 0.094$, and $P_{\delta }\approx 0.409$. Thus, we are able to calculate entropies of the individual cone cells and, as a result, that of $S\left(x_d\right)$ as

$$\begin{array}{cc}\hfill S\left(x_d\right)& =-\sum _i{P}_{i}ln{P}_i\hfill \\ \hfill & \approx 0.324+0.293+0.268+0.222+0.366\hfill \\ \hfill & =1.473.\hfill \end{array}$$

(39)

Now, the entropy difference between differentiated and stem cells reads

$${\langle \Delta S\rangle }_{s\to d}={\langle S\left(x_d\right)-S\left(x_s\right)\rangle }_{s\to d}=-0.136,$$

(40)

where $S\left(x_s\right)$ behaves as the entropy of a uniform distribution, i.e., $S\left(x_s\right)=-{\sum }_{i=1}^5(1/5)ln(1/5)$$=ln5$; this is due to the fact that, for pluripotent stem cells, there are not yet color preferences.

In order for differentiated tissue to maintain its spatial order and integrity, (38) imposes a lower bound on heat dissipation (cell metabolism) as

$${\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}>-{\langle \Delta S\rangle }_{s\to d},$$

(41)

where we have set ${\langle \Delta i\rangle }_{s\to d}\to 0$ for simplicity. Inequality (41) is strongly held for the values obtained here, that is, $9.450\gg 0.136$. This implies that the development of the avian retina is highly robust. In other words, the arrow of time is almost irreversible in this process.

5. Cell Sensing Radius Limits for Robust Tissue Development

In this section, we derive a relationship between total entropy production and cell sensing radius in the particular case of progenitor cells differentiation into the avian cone photoreceptors. We calculate the limits of the sensing radius and, in the parametric space of the LEUP parameters, we suggest the biologically and physically acceptable regions for robust tissue differentiation.

In (32), we have introduced total entropy production as

$$f={\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}+{\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}+{\langle \Delta N\rangle }_{s\to d}-{\langle \Delta i\rangle }_{s\to d}-ln\left[\frac{Z_1}{Z_2}\right],$$

(42)

where thermodynamic properties of system are related to the LEUP quantities. In order to calculate ${\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}$, by assuming the microenvironmental probability distributions as Gaussians, we have

$${\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}=\frac{{\beta }_d}{2}ln\left[2\pi e{\sigma }_d^2\right]-\frac{{\beta }_s}{2}ln\left[2\pi e{\sigma }_s^2\right],$$

(43)

where we have used (29). The above can be simplified to

$${\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}=({\beta }_d-{\beta }_s)ln{\left[2\pi e\right]}^{1/2}+\frac{{\beta }_d}{2}ln{\sigma }_d^2-\frac{{\beta }_s}{2}ln{\sigma }_s^2.$$

(44)

Now, we need to postulate the forms of ${\sigma }_d^2$ and ${\sigma }_s^2$ and their scaling with the corresponding cell sensing radius. Jiao et al. in [61] have found that the differentiated retina mosaic is hyperuniform, that is, ${\sigma }_d^2=V{R}^A$, where $A<D$, and $D=2,3$; on the other hand, for progenitor cells, we assume—(A4)—a Poisson distribution, i.e., ${\sigma }_s^2=U{R}^D$. Here, R is the sensing radius, and V and U are the densities of the cone and progenitor cells neighbors, respectively. Plugging these formulas in (44) leads to

$${\langle \Delta S_{\mathrm{LEUP}}\rangle }_{s\to d}=({\beta }_d-{\beta }_s)ln{\left[2\pi e\right]}^{1/2}+\frac{1}{2}ln\left[\frac{V^{{\beta }_d}}{U^{{\beta }_s}}\right]+\frac{A{\beta }_d-D{\beta }_s}{2}lnR.$$

(45)

Another term of (42) which needs to be dealt with is ${\langle \Delta N\rangle }_{s\to d}$. From (29), we have

$${\langle \Delta N\rangle }_{s\to d}={\langle {\Sigma }_{i=1}^5({\mu }_{i,d}{N}_{i,d}-{\mu }_{i,s}{N}_{i,s})\rangle }_{s\to d},$$

(46)

where i counts different types of cones: green, red, blue, violet, and double. We note that the population at each tissue reads as $N_{i,j}={\rho }_{i,j}{R}^D$, where ${\rho }_{i,j}$ is the density of microenvironment and $j\in \{s,d\}$. Thus, (46) reduces to

$${\langle \Delta N\rangle }_{s\to d}={\langle {\Sigma }_{i=1}^5({\mu }_{i,d}{\rho }_{i,d}-{\mu }_{i,s}{\rho }_{i,s})\rangle }_{s\to d}{R}^D\equiv \Delta \tilde{\mu }{R}^D.$$

(47)

Putting all terms together, total entropy production (42) becomes

$$f(R,{\beta }_d,{\beta }_s)=C_0+C_1({\beta }_d-{\beta }_s)+\frac{1}{2}ln\left[\frac{V^{{\beta }_d}}{U^{{\beta }_s}}\right]+\frac{A{\beta }_d-D{\beta }_s}{2}lnR+\Delta \tilde{\mu }{R}^D,$$

(48)

where we have assumed ${\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}$, $-{\langle \Delta i\rangle }_{s\to d}$, and $-ln(Z_1/Z_2)$ as constants and have grouped them together as $C_0$, and $C_1\equiv (1/2)ln\left(2\pi e\right)$.

Now that we have an explicit formula for total entropy production at our disposal, we can obtain the optimal sensing radius for which total entropy production reaches its extremum. To this end, we note that the first derivative of (48) vanishes at

$$R_c={\left(\frac{-\Gamma }{\Delta \tilde{\mu }D}\right)}^{1/D},$$

(49)

where $\Gamma \equiv (A{\beta }_d-D{\beta }_s)/2$. When $R_c$ is a positive quantity, then the inequality $\Gamma /\Delta \tilde{\mu }<0$ should always be satisfied. In order to determine the conditions for which $R_c$ minimizes or maximizes total entropy production, we calculate the second derivative of (48) at (49), and obtain

$$\frac{{\partial }^{2}f}{\partial R^2}{|}_{R_c}=-\frac{\Gamma D}{R_c^2}.$$

(50)

Thus, $R_c$ minimizes total entropy production if $\Gamma <0$ and it maximizes if $\Gamma >0$.

In the left panel of Figure 5, we have illustrated total entropy production (48) as a function of R for a particular set of the LEUP parameters which minimize total entropy production based on (50). The values of the curve above the R-axis (gray line) are biologically and physically relevant as total entropy production is positive, i.e., it ensures the robustness of the differentiation process. Moreover, this plot shows that, as we move away from $R_c$, (positive) total entropy production is rapidly increasing. The right panel of the figure illustrates f as a function of R and ${\beta }_d$, where the acceptable regions lie above the gray plane. Interestingly, minimization of total entropy production occurs only for $\Delta \tilde{\mu }>0$, which assumes a decrease in the proliferative activity of differentiated phenotypes.

In Figure 6, the case for which total entropy production reaches its maximum is shown. This figure illustrates that (positive) total entropy production is bounded in this case. The immediate implication of this is that tissue dedifferentiation is possible (negative total entropy production) for a large sensing radius. If we want to have a positive total entropy production, in order to avoid reversibility of differentiated tissue, we must have some sort of fine-tuning in order to restrict the values of R in such a way that they lead to a positive f. According to Bialek’s postulated biophysical principles [30], fine-tuning in Nature is not favorable. In the case of the avian photoreceptors, we can further constrain our parameters and identify further arguments to exclude this case.

Application: The Average Sensing Radius of the Avian Cone Cell

In this subsection, we calculate the average sensing radius for each photoreceptor, and, by exploiting (48), we look for the regions, in the parametric space of $({\beta }_s,{\beta }_d)$ and $({\beta }_s,{\beta }_d,\Delta \tilde{\mu })$, which result in a positive total entropy production and ensure the robustness of the differentiated tissue.

In the retina of a bird, as explained in the previous sections, the five types of photoreceptor cells, namely, green, red, blue, violet, and double cones, form mosaic structures. From the experimental data (see Supporting Information of [60]), we have information about the average standard deviation of the Nearest Neighbor Distribution (NND) for each color as: ${\sigma }_g\approx 1.248$, ${\sigma }_r\approx 1.548$, ${\sigma }_b\approx 1.729$, ${\sigma }_v\approx 2.292$, and ${\sigma }_{\delta }\approx 0.948$. In addition, Jiao et al. [61] have found how the variance of an avian cone ${\sigma }^2$ is related to its sensing radius R as

$${\sigma }^2\left(R\right)=M_1{R}^2+M_{2}RlnR+M_{3}R,$$

(51)

where $M_1$, $M_2$, and $M_3$ have been calculated for each color (see Table I of Ref. [61]). By having this information at our disposal, we are able to calculate the average sensing radius for each photoreceptor as

$${\overline{R}}_{\delta }\approx 0.786,{\overline{R}}_g\approx 1.426,{\overline{R}}_r\approx 2.006,{\overline{R}}_b\approx 2.288,{\overline{R}}_v\approx 3.553.$$

(52)

We note that, if we assume that the minimum sensing radius is equal to the cone size, then the above sensing radii are very close to the average cone size $2.59\pm 1.05\mathsf{\mu }$m (here, we have used the $\pm 3\sigma$ rule for the data presented in Figure 5A of Kram et al. [60]).

In the beginning of this section, we have introduced the variance of the differentiated cells as ${\sigma }_d^2=V{R}^A$; now, by approximating (51) to have such a particular form and setting $V=1$, and using the values of (52), we can find A for each cone as

$$A_{\delta }\approx 0.440,A_g\approx 1.247,A_r\approx 1.256,A_b\approx 1.322,A_v\approx 1.308,$$

(53)

where they are in agreement with the assumption of hyperuniformity made previously, that is, $A<2$.

As an illustration, based on the values of (52) and (53), we have demonstrated the biologically (and physically) acceptable regions of total entropy production (48) in the parametric space of the LEUP parameters: $({\beta }_s,{\beta }_d)$ and $({\beta }_s,{\beta }_d,\Delta \tilde{\mu })$ in Figure 7, in the case of the double cone photoreceptor.

At this point, we should remind that the parameter $\Delta \tilde{\mu }$ is proportional to the difference of differentiated and pluripotent cells division times. It is well known that the division time of any differentiated cell is much larger than the division time of the corresponding progenitor one [64]. Therefore, it is safe to assume that $\Delta \tilde{\mu }>0$. This allows us to identify a relationship for the corresponding $\beta$’s of the cones and their progenitors, where Figure 7 implies that ${\beta }_d>{\beta }_s$. The latter is a sensible result since differentiated cells are expected to be more attentive to their microenvironment in order to ensure optimal cooperation and tissue integrity.

6. Discussion

In this paper, we have posed the question (Q1) on how cells coordinate intrinsic and extrinsic variables to determine cell decisions that eventually lead to organized and stable tissues. To tackle this problem, we have employed the Least microEnvironmental Uncertainty Principle (LEUP), which has been recently proposed to understand cell decision-making in multicellular systems, and so far it has been applied to cell migration force distribution [25], collective cell migration [26], and binary phenotypic plasticity [27]. In the context of the LEUP, we regard differentiation as a sort of Bayesian decision-making, where cells update their intrinsic variables by encoding microenvironmental information and producing relevant responses. This provides us with a distribution of internal states that depends explicitly on the information of the cell current microenvironment, which is represented by a mesoscopic microenvironmental entropy. Interestingly, we have shown that local microenvironmental entropy should decrease in time leading to more organized cellular microenvironment, which is the case in differentiated tissues. As a proof of principle, we have challenged the LEUP predictions to reproduce differentiated avian photoreceptor mosaics. Although, by fitting a single parameter $\beta$, we have successfully reproduced the photoreceptor statistics, this still cannot be considered as a rigorous validation. To this end, we have recently gathered an inter-species collection of photoreceptor mosaics to further investigate the potential of the LEUP to reproduce these tissues and possibly classify them.

By using the aforementioned results, we have attempted to shed light on the macroscopic transition between pluripotent and differentiated tissues and have specified it to the formation of photoreceptor mosaics, which is related to the question (Q2) posed in the Introduction. In this respect, we have developed a stochastic thermodynamic-like theory, based on the Crooks’ theorem, for a general cell and tissue differentiation process. We have shown that differentiated tissues are highly robust to dedifferentiation, even though individual cells are allowed to go back into pluripotent phenotypes. Biologically, the robustness of differentiated tissues depends on reduced proliferation, a change from anaerobic to aerobic metabolism, and increased cell sensing that leads to a higher order of microenvironmental organization. In particular, we have estimated the critical sensing radii of photoreceptor cones that ensure the thermodynamic robustness of differentiated mosaics, which turns out to be in the range of $0.8$ $\mathsf{\mu }$m to $3.5$ $\mathsf{\mu }$m, see (52). We note that the critical radius is the minimal radius required to ensure tissue robustness and therefore should serve as a lower bound for the real values. Now, if we assume that the minimum sensing radius is equal to the cone size, then our predicted range correlates with the average cone size $2.59\pm 1.05\mathsf{\mu }$m (here, we have used the $\pm 3\sigma$ rule for the data presented in Figure 5A of Kram et al. [60]).

In summary, our LEUP-driven model is based on the four crucial assumptions: (A1) there is a timescale separation between the internal and microenvironmental variables’ dynamics, (A2) the multicellular system (tissue), where the cell is differentiating, follows a Markovian dynamics with the assumption of microscopic reversibility, (A3) a flat cell state distribution is assumed for pluripotent cell states, and (A4) the spatial distribution of the early microenvironmental pluripotent cells follows a Poisson distribution. Based on these assumptions, we have arrived at three important results: (I) predicting the color percentage of the cone cells in the avian retina without any knowledge about the underlying biophysical and biochemical mechanisms, (II) demonstrating the robustness of cell-tissue differentiation in thermodynamic terms, and (III) determining the limits of the cell sensing radius by establishing a relation between total entropy production and microenvironmental sensing. In the following, we further elaborate on these results.

(I) Prediction of the cone color distribution: By calibrating a single parameter in the LEUP, we are able to predict the cone color percentage in the avian retina accurately. Our finding regarding the LEUP parameter, which reads as $\beta \approx 1.754$ (close to 2), gives a strong indication that cells sense their environment quasi-optimally when choosing a particular cell fate during the differentiation process. For future study and investigation, we want to examine the validity of this result for photoreceptor mosaics of other species [65,66,67,68].

(II) Robustness and measurement of information gain in differentiation: We have constructed the fluctuation theorem for tissue differentiation and have derived a generalization of the second law of thermodynamics for this process based on Markovian dynamics. It would be interesting and more realistic to relax this assumption and to analyze the problem using the LEUP on the basis of non-Markovian or memory processes [69], which can lead to different results. The only requirement is that system should have a unique stationary state. We also note that the LEUP-driven second law of thermodynamics can also be seen as a generalization of the Bayesian second law of thermodynamics [70] and the conditional second law of thermodynamics in a strongly coupled system [71].

Our theory, which in the current paper has been applied to the specific case of the avian photoreceptor mosaics, suggests that differentiated tissue is (highly) thermodynamically robust, that is, the arrow of time is almost irreversible, and this robustness depends on microenvironmental sensing and cell metabolism. It should be remarked that, as is demonstrated in (38), if we have the values of ${\langle {\beta }^{{}^{\prime }}\Delta Q_{x_s\to x_d}^{\tau }\rangle }_{s\to d}$ and ${\langle \Delta S\rangle }_{s\to d}$, then we can determine the upper bound of pointwise mutual information difference which is denoted as ${\langle \Delta i\rangle }_{s\to d}$ within the context of the LEUP. However, in the present work, we have only considered transitions between equilibrium end-states and have set ${\langle \Delta i\rangle }_{s\to d}=0$. In a future study, we plan to investigate non-equilibrium dynamics of transitions between progenitor and differentiated cell states and establish the upper bound of information gain in cell differentiation. (We also note that the value of ${\langle \Delta i\rangle }_{s\to d}$ can also be obtained directly from experiments, see Refs. [72,73,74,75].)

(III) Limits of sensing radius: By studying total entropy production as a function of cell sensing radius and the LEUP parameters, we have provided an understanding of how a cell regulates its sensing radius according to its microenvironment to ensure the thermodynamic robustness of differentiated tissue. We have shown two cases where (a) the entropy production goes to infinity beyond a certain threshold radius which is depicted in Figure 5 and (b) the entropy production goes to a maximum value as in Figure 6. Interestingly, we conclude that the former is the most biologically relevant case since it requires the division time of differentiated cells to be larger than that of the pluripotent ones and biologically systems operate away from fine-tuned parameter regimes to withstand noisy perturbations. On the technical side, we have assumed that the spatial distribution of pluripotent tissue resembles a Poisson distribution. It would be interesting to relax this assumption and to derive this distribution from real tissue data.

One important issue is the range of validity of (A1) regarding the timescale separation between cell decision and cell cycle characteristic times. Although cell decisions may seem to be happening within one cell cycle, the underlying molecular expressions may evolve over many cell cycles [76,77]. When these molecular expressions cross a threshold, then cell decision emerges very fast. Therefore, the definition of cell decision should be treated with care. In our case, we specify a cell decision only when the cell state switches to another dynamic attractor that induces at the same time some noticeable phenotypic changes. Such attractor transitions are manifested as switches with much shorter characteristic times than a cell cycle [20].

We want to make a brief comment on the relation of our theory to the commonly used approach of the maximum entropy production (MEP). The MEP formalism explains only the transition from a pluripotent state to a differentiating one without realizing the corresponding dynamics. To build a connection between these two, one has to construct the LEUP theory for transition paths like maximum caliber principle [78]. Instead of internal and external variables, one uses internal and external paths of the corresponding evolution. Then, by exploiting the formulation of maximum caliber, one can write the time evolution of microenvironmental path entropy—as a conservation equation—in terms of sources and fluxes and subsequently in terms of path action and entropy production [79,80]. In this regard, one can use maximum caliber principle to construct appropriate transition probabilities and even understand the spatiotemporal dynamics [81]. Finally, one should maximize the internal path entropy which resembles the MEP approach. Working out the details of this connection remains for a future work.

Our proposed theory has important and interesting implications for cancer research and therapy. In particular, (42) states that the balance of metabolic, proliferative, and tissue organization changes (the LEUP term), needing to take place in order to destabilize the differentiated state, that is, to promote carcinogenesis. Until now, the majority of the therapies were focused on antiproliferative strategies, such as chemotherapy and radiotherapy, and more seldom to the metabolic conditions such as vasculature normalization. Here, we have proposed that changes in the tissue organization play a critical role. This fact has been very recently identified in the context of tumor evolution by West et al. [82]. We also note that, in [83], it is realized that microenvironment normalization might be the key for immunotherapeutic success. The mechanistic connection between tissue architecture and cell sensing mechanisms is established in the context of our theory. Strikingly, the experimental work of Levin’s group [84] shows that disrupting the ion channel sensing in a tissue can induce tumorigenesis. In this regard, we have put forward that investigating changes in the cell sensory processes deserve more attention and might be pivotal in treating cancers. Our goal is to calibrate the existing theory to human photoreceptors data, thus we could apply these ideas to retinoblastoma tumours.

In a nutshell, we have shown how the LEUP facilitates the inference of cellular intrinsic states (or, cell phenotypes) by means of local microenvironmental entropies or fluctuations. This allows the evaluation of cellular states without the detailed grasp of the underlying mechanisms. The sole knowledge about extrinsic variables distributions (or, collective cell decision-making) suffices. Therefore, we can apply the LEUP to cell differentiation problems where the biological or biophysical knowledge is unclear or unknown.