
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01794-2
10.1016/j.isci.2024.110569
110569
Article
Microenvironmental entropy dynamics analysis reveals novel insights into Notch-Delta-Jagged decision-making mechanism
Pujar Aditi Ajith 1269
Barua Arnab 39
Dey Partha Sarathi 1
Singh Divyoj 127
Roy Ushasi 18
Jolly Mohit Kumar mkjolly@iisc.ac.in
1∗
Hatzikirou Haralampos haralampos.hatzikirou@ku.ac.ae
4510∗∗
1 Department of Bioengineering, Indian Institute of Science, Bangalore 560012, India
2 Undergraduate Program, Indian Institute of Science, Bangalore 560012, India
3 Tata Institute of Fundamental Research, Hyderabad 500046, India
4 Mathematics Department, Khalifa University, P.O. Box: 127788, Abu Dhabi, UAE
5 Technische Univesität Dresden, Center for Information Services and High Performance Computing, Nöthnitzer Straße 46, P.O. Box: 01062, Dresden, Germany
∗ Corresponding author mkjolly@iisc.ac.in
∗∗ Corresponding author haralampos.hatzikirou@ku.ac.ae
6 Present address: Department of Physics, The University of Texas at Austin, Austin, TX 78712, USA

7 Present address: Department of Physics, University of California, Santa Barbara, CA 93106, USA

8 Present address: Department of Physics, Indian Institute of Science Education and Research, Pune 411008, India

9 These authors contributed equally

10 Lead contact

21 8 2024
20 9 2024
21 8 2024
27 9 11056919 3 2024
31 5 2024
19 7 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by/4.0/ This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Summary

Notch-Delta-Jagged (NDJ) signaling among neighboring cells contributes crucially to spatiotemporal pattern formation and developmental decision-making. Despite numerous detailed mathematical models, their high-dimensionality parametric space limits analytical treatment, especially regarding local microenvironmental fluctuations. Using the low-dimensional dynamics of the recently postulated least microenvironmental uncertainty principle (LEUP) framework, we showcase how the LEUP formalism recapitulates a noisy NDJ spatial patterning. Our LEUP simulations show that local phenotypic entropy increases for lateral inhibition but decreases for lateral induction. This distinction allows us to identify a critical parameter that captures the transition from a Notch-Delta-driven lateral inhibition to a Notch-Jagged-driven lateral induction phenomenon and suggests random phenotypic patterning in the case of lack of dominance of either Notch-Delta or Notch-Jagged signaling. Our results enable an analytical treatment to map the high-dimensional dynamics of NDJ signaling on tissue-level patterning and can possibly be generalized to decode operating principles of collective cellular decision-making.

Graphical abstract

Highlights

• LEUP-based low-dimensional modeling recapitulates NDJ spatiotemporal patterns

• Abstraction of the NDJ pathway into a cell decision-making process

• Identifying a key parameter for transitioning from lateral inhibition to induction

• Analytical results map mechanistic pathway models to the LEUP framework

Molecular network; Biocomputational method; Flux data; In silico biology; Biological constraints

Subject areas

Molecular network
Biocomputational method
Flux data
In silico biology
Biological constraints
Published: August 21, 2024
==== Body
pmcIntroduction

The Notch signaling pathway is evolutionarily highly conserved across species and is crucial during many steps of embryonic development as well as in maintaining tissue homeostasis during adulthood.1 It is a cell-cell communication pathway and contributes to cellular proliferation, differentiation, and apoptosis.2,3 Thus, its dysregulation leads to various developmental defects and can also be critical in cancer progression.1 In cancer, it can promote metastasis and the associated phenomenon of epithelial-mesenchymal transition (EMT).1,4

The Notch signaling pathway is juxtacrine, i.e., the extracellular domain of the Notch receptor of one cell binds to Delta and/or Jagged family ligands of the neighboring cells (trans-activation). This binding enables the proteolytic cleavage of Notch, whose intracellular domain—Notch intracellular domain (NICD)—then trans-locates to the nucleus to initiate downstream signaling. Notch-Delta signaling between two neighboring cells leads to lateral inhibition, where two cells adopt opposite fates—one acts as a sender (high Delta, low Notch) and the other as receiver (high Notch, low Delta).5 However, Notch-Jagged signaling between two cells often leads to lateral induction, where both cells acquire a similar phenotype of being a hybrid sender/receiver (high Notch, high Jagged). Thus, the spatiotemporal dynamics of Notch-Delta-Jagged signaling contributes to various instances of tissue-level pattern formation and has attracted significant attention in computational modeling.6,7,8

Most computational models developed for Notch signaling so far have incorporated multiple parameters and are high dimensional in nature in terms of variables (Delta, Notch, Jagged, NICD, glycosyltransferase Fringe, NICD targets etc.) considered,8,9,10,11,12,13,14 thus making it difficult for any analytical treatment of the same. Thus, identifying the fundamental design principles of Notch signaling-driven cellular decision-making in a low-dimensional space remains challenging.

Cellular decision-making incorporates various factors: (1) internal components of the cell (metabolite levels, structure of gene regulatory network, etc.), (2) external components (biochemical and biophysical microenvironment of the cell, density of neighboring cells, etc.), and (3) noise or uncertainty in sensing the cellular microenvironment, and in various biochemical processes happening inside the cell.15 Although single-cell decision-making has been studied meticulously, the study of cellular decision-making on a multi-cellular level remains obscure. Although from a different approach, the intracellular and inter-cellular communications have been understood using the Hopefield network.16,17,18 To address this, we recently proposed a principle known as the least micro-environmental uncertainty principle (LEUP).19 The principle is based on the statistical mechanical premise, which can be applied in different biological contexts: collective cell migration, cell differentiation, cellular sensing dynamics, and go-or-grow dynamics.20,21,22,23 It adopts the idea of Bayesian inference to construct a mathematical framework. Although there are high-throughput technologies that give us static pictures of a cell’s internal variables, the knowledge of the dynamics of those variables is partial.23 LEUP attempts to compensate for this missing knowledge by employing a variational principle on the probability distribution of a cell’s microenvironment which arises from two factors: (1) a cell’s uncertainty in sensing its microenvironment and (2) the phenotypic distribution in the cell population. This is codified in a microenvironmental entropy term. LEUP asserts that the variational principle leads to an ordinary differential equation (ODE) description of the dynamics of the internal cellular state in terms of microenvironmental entropy. The inspiration behind this principle comes from Friston’s free energy principle24 and Bayesian brain hypothesis.25 Analogous concepts have also been proposed in the seminal works of William Bialek.26

In this work, we introduce a sensitivity parameter that quantifies a cell’s capacity for sensing its environment (through ligand binding, mechano-sensing, pseudopodia extension, ion channels, etc.). We derive the relationship between this sensitivity parameter and the dynamics of the cell’s internal and external variables. By demonstrating how cellular behavior depends on the ratio of the biophysical forces and variance of the microenvironmental probability distribution, we predict the phenotypic distribution of a population of cells from varying degrees of information about cell-cell interaction and cell-autonomous variables.

In the context of the Notch signaling pathway, we show that the tissue-level patterning due to Delta-driven lateral inhibition and Jagged-driven lateral induction can be described in terms of the time evolution of local entropy in different regimes of the sensitivity parameter. Specifically, in lateral inhibition (induction), the sensitivity parameter has a negative (positive) value and the local entropy increases (decreases) over time. In scenarios where neither Notch-Delta nor Notch-Jagged signaling dominates, the parameter approaches zero, resulting in the emergence of random patterns.

Hence, we have been able to map the high-dimensional mathematical model to a low-dimensional description amenable to analytical treatment. Our formalism can be generalized for elucidating the operating principles governing collective cellular decision-making processes. Additionally, it holds promise for providing novel insights into specific biological phenomena, as exemplified in this study.

This article has been arranged in the following way: in STAR Methods, we review the Bayesian decision-making and the derivation of the phenotypic Langevin equation. Furthermore, we calculate the sensing intensity parameter using a small-noise approximation applied to the Notch-Delta-Jagged model. After that, in Section results we study the pattern formation and phase transition regimes for different values of sensitivity β and interaction radius r in square lattices and later compare with the mechanistic model of the Notch-Delta-Jagged system. Finally, we discuss and conclude the outcomes of our analyses in Section discussion.

Results

At first, we shall summarize the local microenvironmental information for the Notch-Delta-Jagged system, and then we shall compare the behavior of the LEUP-driven Langevin equation on the lattice. The schematics have been described in Figure 1.Figure 1 Schematic diagrams representing the two systems studied here

(A) This illustrates the Notch-Delta-Jagged signaling pathway. The trans-membrane protein Notch (N) can undergo trans-activation with ligands Delta (D) and Jagged (J) releasing NICD (I), a protein marker that then inhibits D expression and activates N and J production. When both components of this receptor protein-ligand complex are from the same cell, they simply degrade without NICD release, i.e., they undergo cis-inhibition. This schematic is inspired by Boareto et al.27

(B) This represents the LEUP-based phenomenological model which has only two parameters—a radius of interaction r and sensitivity β. Each cell senses its internal state (phenotype), Xn, and external variables (Yn’s) which are the phenotypes of its neighboring cells.

(C) The principles governing LEUP are illustrated. S(X), the entropy of the cell population’s phenotype, reduces over time, which can be visualized as phase space volume over X1,…,XN decreasing over time.

Sensitivity β and interaction radius r induce phase transitions with respect to microenvironmental entropy

The first question that we attempt to answer is how sensitivity parameter β and interaction radius r impact the dynamics of microenvironmental entropy and pattern formation in LEUP spatiotemporal simulations. Despite the low dimensionality of the parametric space, only r and β, the LEUP systems display a wealth of behaviors.

First, we investigate the impact of the parameters on the steady-state distribution of the phenotypic parameter Xnt. We obtain these distributions by constructing the histogram of phenotypes over all lattice nodes. For negative sensitivity β and nearest neighborhood interactions r=1, we observe the emergence of multimodal distributions. In particular, for decreasing β<0 a transition from bimodality to trimodality occurs. Interestingly, this effect is lost for higher radii such as r=5. In the β>0 case, we observe a highly concentrated distribution around the Xnt=0 value. This implies that the majority of lattice nodes acquire a “hybrid” phenotype. These results can be readily seen in Figure 2A.Figure 2 Phenotype distribution and the average rate of change of microenvironment uncertainty over time are examined

(A) Kernel density estimates (histograms smoothened with Gaussian kernels) are plotted for different LEUP cell populations. The x axis is the range of phenotypes, Xnt∈[−1,1], and the y axis is the phenotype distribution, mostly normalized to be in [0,1], exceeding 1 in cases where the peaks are too sharp (β>0 populations). As |β| increases (i) (r=1,β<0) populations go from bimodal to trimodal, showing increased sensitivity selects for the hybrid Xnt=0 state. In (ii) and (iv), the existence of a threshold sensitivity |βthresh| is clearly delineated, only beyond which cell populations are unimodal/very sharply peaked instead of uniformly randomly distributed.

(B) These plot the average slope of ⟨σ2⟩L over time, i.e., ⟨d⟨σ2⟩Ldt⟩T averaged over time, for each β for different r. For (i) r=1, βthresh− is observed and marked by (a). It is the sensitivity beyond which a noticeable increase in ⟨d⟨σ2⟩Ldt⟩T shows the cells beginning to sense. (ii) r=5. Here, both βthresh− and βthresh+ are marked by a decrease in ⟨d⟨σ2⟩Ldt⟩T, denoted by (a) and (b). We can see that the transition to informed decision-making is more sharp for β>0, evidenced by the sharp drop in ⟨d⟨σ2⟩Ldt⟩T. Throughout this figure, blue (green) plots correspond to negative (positive) β.

Now, we focus on the impact of β on the local variance growth rate, which is a proxy of a temporal average of microenvironmental entropy rate, for difference sensing radii (Figure 2B). When cells sense the immediate neighborhood r=1, for negative sensitivity we observe positive growth rates of the observable ⟨d⟨σ2⟩Ldt⟩T. This is consistent with the multi-stable steady-state distributions. On the contrary, for β>0 the local variance decays in time implying that the steady state approximates a delta distribution. The latter is consistent for larger radii for instance for r=5. However, for large sensing radii the ⟨d⟨σ2⟩Ldt⟩T stays non-positive for variations of β. Also, the behavior is more complex exhibiting a non-monotonic behavior.

To further investigate the latter observation, we plot the steady-state local variance ⟨σ2⟩L for a large combination of sensitivity parameter β and the radius r in Figure 3 (also see Figure S1). Please note that the local variance of microenvironmental phenotypes has been calculated using the Moore neighborhood over all the lattice nodes and the averages have been taken over the ensembles of different realization close to the steady states only. Details of these simulations can be found in Figure S1. For any radius and for large enough β>βthresh+, we observe always an unimodal delta distribution of minimal entropy. For radii r>1 and sensitivities βthresh−≤β≤βthresh+, the local variance stays high. However, when β crosses the lower bound βthresh−, we observe again a microenvironmental variance reduction. Interestingly, for r=1 case, there is simply an abrupt phase transition from high to low ⟨σ2⟩L when β crosses threshold βthresh+.Figure 3 The phase space diagram of the average microenvironmnental uncertainty as a function of sensing radius and sensitivity

The steady-state ⟨σ2⟩L values are plotted across the phase space of the LEUP model. The increase in |βthresh| as r increases is very evident, marked by (i) for β>0 and (ii) for β<0. A third gray arrow (iii) shows the increase (corresponding to the cells beginning to sense for sensitivity beyond βthresh−(r=1)) and gradual decrease in ⟨σ2⟩L as trimodality emerges.

Lateral inhibition (induction) leads to increasing (decreasing) microenvironmental entropy dynamics

The Notch-Delta-Jagged signaling system has been extensively studied and analyzed.4,27 However, the NDJ system has not been analyzed on the basis of local entropic dynamics. Here, we conduct spatial simulations of Equation 17 to understand the behavior of the microenvironmental entropy, in terms of local variance ⟨σ2⟩L (as introduced in the STAR Methods). We consider the NICD levels to correspond to the phenotype of the cells. In turn, we plotted ⟨σ2⟩L as a function of time for the Delta-dominated lateral inhibition regime (Figure 4A(i)) and the Jagged-dominated lateral induction regime (Figure 4B(i)).Figure 4 The two Notch-Delta-Jagged signaling regimes are shown

(A) Lateral inhibition (checkerboard spatial patterning) and (B) lateral induction (uniform spatial patterning) are observed in (i) Notch-Delta heatmaps of NICD levels. These are normalized to be between [0,1]. Besides patterning, the time evolution of ⟨σ2⟩L has also been plotted for both cases. It increases for lateral inhibition (A(ii)) and decreases for lateral induction (B(ii)) over time.

The simulations have shown that, for lateral inhibition, ⟨σ2⟩L increased over time as shown in Figure 4A(ii). For the lateral induction regime, where Jagged is overexpressed, we observe a decrease in the local variance over time (Figure 4B(ii)). This indicated that the NDJ system agrees with the LEUP premise that microenvironmental should decrease or increase for certain parameter constellations. In the following, we delve deeper into the connection between local entropy and spatial dynamics.

LEUP recapitulates a noisy NDJ spatial patterning

The Notch-Delta-Jagged mechanism induces distinct patterns when Jagged is overexpressed or not, as seen in Figure 4. In particular, the lateral induction dynamics result into a homogeneous spatial distribution of hybrid phenotype, whereas the lateral inhibition induces the infamous checkerboard-like patterns. The question is if an LEUP model, with only two parameters, is able to capture these two equilibrium patterns.

Regarding the lateral induction case, we observe, when β>βthresh+, a uniform distribution of phenotypes Xn=0 emerges. This corresponds to a uniform pattern exactly like the one of the Jagged-dominated NDJ system. In the case of lateral inhibition, things are not so trivial. The NDJ system typically produces frustrated chessboard patterns, whereas the LEUP simulation produces a noisy version of them. In order to quantitatively compare NDJ and LEUP patterns, we invoke the radial distribution function (rdf) and the spatial fast Fourier transform (FFT). In Figure 5C, we observe that the NDJ is close to the ideal chessboard rdf and the LEUP for r=1 and β<0 is a noisy version of that. The same conclusion can be drawn by observing Figure 5B(i and iii), where the FFT and the corresponding power spectrum show the LEUP and NDJ chessboard-like versions.Figure 5 The phenotypic map of the cells is studied in terms of radial distribution function

(A) Heatmap of LEUP phenotypes for (r=2,β=−100) shows isolated “islands” of extrema cells (Xnt=±1) surrounded by hybrid cells (Xnt=0).

(B–D) (B) Heatmaps of spatial Fourier transforms of (i) chessboard patterned LEUP (r=1,β=−100), (ii) LEUP (r=2,β=−100), and (iii) chessboard-patterned NDJ (D0,J0)=(1600,700) are plotted. (i) and (iii) correspond to similar spatial patterns and are also alike, with high contributions from a few low frequencies due to spatial ordering. However, the ring-like structure present in (ii) shows that non-trivial spatial order exists in the (r>1,β<0) systems as well. To characterize the spatial patterning, the (C) radial distribution functions and the (D) radial power spectra of the three lattice systems are plotted, along with those of a perfect chessboard pattern for reference. While (i) and (iii) showed reasonably good correspondence with the reference chessboard as expected, the structure in (ii) shows in the power spectra. Whereas the remaining power spectra peaked at low values of r in the frequency domain, (corresponding to low frequencies), LEUP (r=2,β=−100) systems peak at a noticeably larger radius (in frequency domain), corresponding to larger frequencies.

Now let us focus on the patterns produced by LEUP when the sensing area has radius r>1. Interestingly, the chessboard pattern is not any more attainable. However, we are obtaining a sparse activation of nodes. In particular, cells expressing extreme phenotypes are almost always surrounded by cells of hybrid (Xnt=0) phenotypes; rarely, there are more than a couple of Xnt=±1 cells as immediate neighbors. The FFT and the rdf indicate that phenotype activation takes place further away from the immediate neighborhood (Figures 5B(ii), 5C, and 5D). These patterns are called “checkerboard” in the literature. Interestingly, when the NDJ system communicates in larger radii, then the patterns are similar to those observed in developing pancreas patterning, in which the Notch-Delta signaling adopts a lateral stabilization28 mechanism.

The sensitivity parameter (β) quantifies the relative impact of intrinsic and extrinsic variables on cell decisions

The sensitivity β quantifies a cell’s ability to sense and integrate information about its microenvironment. Now, we focus on the question of calculating β when the sensing biophysical processes are known and well described by a set of differential equations. To answer this question, for the n− th cell, we consider the external variables Yns, such as ligand densities, and internal variables, Xns, such as receptors or internalized proteins. Let us assume a dynamical system of internal and external variables as(Equation 1) X.nt=G(Xnt,Ynt)+ξn,Y.nt=F(Xnt,Ynt)+ηn,

where the intrinsic flow can be rewritten as G(Xnt,Ynt)=G˜(Xnt,Ynt)−Xnt and the extrinsic one as F(Xnt,Ynt)=F˜(Xnt,Ynt)−Xnt. To simplify the calculations, we assume the steady-state responses of internal and external variables as(Equation 2) Xns=G˜(Xns,Yns)+ξn,Yns=F˜(Xns,Yns)+ηn,

where the Xns and Yns are the corresponding steady states of intrinsic and extrinsic variables, respectively. Here, we assume that the aforementioned noises follow the following multivariate Gaussian distributions ξn∼N(0,ΣXns) and ηn∼N(0,ΣYns). The probability distribution of internal variables according to the aforementioned dynamics reads as(Equation 3) P(Xns∣Yns)=e−12(Xns−G˜(Xns,Yns))ΣXns−1(Xns−G˜(Xns,Yns))T(2π)d2∣ΣXns∣12.

In the small-noise approximation, using similar arguments with Bialek,26 we can calculate the connection between the internal and external variable probability distributions (akin to a change of variables during integration).(Equation 4) P(Xns)||dXnsdYns||=P(Yns).

The quantity ||dXnsdYns|| corresponds to the absolute value of the determinant of the Jacobian |∇YnsG(Xns,Yns)| estimated at the steady state (please note that ∇YnsG(Xns,Yns)=∇YnsG˜(Xns,Yns)). We substitute the value of Jacobian in the small-noise approximation limit and the value of P(Yns) into Equation 4:(Equation 5) P(Xns)|∇YG(Xns,Yns)|=e−12(Yns−F˜(Xns,Yns))ΣYns−1(Yns−F˜(Xns,Yns))T(2π)d2∣ΣYns∣12

The quadratic term related to the external dynamics Gaussian goes to zero because we expect that the random variable Yns is very close to its expected value F˜(Xns,Yns) at steady state. At this point, we use the LEUP steady state as calculated in Equation 12. Entropy of the microenvironment, a d-dimensional multivariate normal distribution, is 12ln((2πe)d∣ΣYns∣Xns∣). Thus,(Equation 6) P(Xns)=e−β2ln((2πe)d∣ΣYns∣Xns∣)Z=∣ΣYns∣Xns∣−β/2Z′

where the normalization constants are Z=∫…∫e−βS(Yns∣Xns)dXns and Z′=∫…∫∣ΣYns∣Xns∣−β/2dXns. Substituting into Equation 5, we get(Equation 7) |ΣYns∣Xns|−β/2Z′=|∇YG(Xns,Yns)|(2π)d2∣ΣYns∣12

Taking the logarithm and simplifying gives us(Equation 8) −β2ln|ΣYns∣Xns|=−ln|∇YG(Xns,Yns)|−ln((2π)d/2∣ΣYns∣1/2Z′)

Here, we estimate the normalization term as Z′=|ΣYns|−β/2M, where M=∫dXns is the capacity concentration of intrinsic variables. Also, we use the definition of entropy to obtain ln((2π)d/2|ΣYns|1/2)=S(Yns)−d/2. We further simplify this equation (see sensitivity parameter β calculation for details), and we can finally arrive to write sensitivity or β as(Equation 9) β=−ln|∇YG(Xns,Yns)|+S(Yns)−cI(Xns:Yns)

where c=lnM+d/2 is a constant scalar and the mutual information for multivariate Gaussian reads I(Xns:Yns)=12ln(|ΣYns||ΣYns∣Xns|).

Sensitivity β encodes transitions from uni- to bimodality in Notch-Delta-Jagged phenotypic distributions

In the previous section, we have derived an analytical relation for the sensitivity parameter β. In particular, it informs us about the phenotypic responses on microenvironmental variations scaled by the degree of intrinsic-extrinsic correlations. Here, we investigate the relevance of β in the context of NDJ phenotypic dynamics.

As we have seen before, the sensitivity parameter controls the transitions from unimodality to multimodality of phenotypic probability distribution function in LEUP simulations. In particular, when β<0, the phenotypic distribution is bimodal or else unimodal. In turn, we simulate the spatial Notch-Delta-Jagged system for varying (D0,J0) while keeping all other parameters constant. These parameters represent the basis production rates of Delta and Jagged ligands and are the primary drivers of either lateral inhibition (Delta-dominated) or lateral induction (Jagged-dominated) cell signaling (as shown by Boareto et al.4). We ran simulations for plausible ranges of these production rate, i.e., (D0,J0)∈[500,2000]2, subsequently calculated β for each simulation, from Equation 8, where details can be found in the STAR Methods section, and calculated the corresponding modality index (as defined in the STAR Methods section).

In Figure 6, there is a great correspondence between the phase diagrams of β and the modality index. The simulations that belong to the (D0,J0) bimodality regime correspond to a negative β, and the rest of the phase space indicates unimodality and β<0. The parameter β even captures some non-linear features of the landscape, and the phase space boundary between induction and inhibition bears a remarkably close resemblance between the two. This is more evident in Figure 7A where the boundary of the regimes can be clearly indicated. Finally, we can view β as a quantitative mapping between the LEUP cell decision-making model and the mechanistic NDJ signaling model.Figure 6 Visualizing sensitivity parameter and modality index in (D0,J0) phase space

Here we plot the values of (i) β and (ii) modality index, defined in Equation 22 across the (D0,J0) phase space of the Notch-Delta-Jagged signaling pathway. Negative (positive) values of both order parameters corresponding to bimodal (unimodal) populations are highlighted. While the former is a calculated quantity capturing features of the phase space, the latter is a statistical index that ratifies it.

Figure 7 Effect of dynamics of extrinsic and intrinsic variables on phenotypic response

(A) The phase diagrams of (i) β and (ii) modality index are discretized by naive thresholding, i.e., all positive (negative) values are equated to +1(−1). This illustrates the phase boundary between the positive and negative regimes.

(B) The interpolated phase diagram of a quantity related to the dynamics, log|∇XG(Xns,Yns)∇YG(Xns,Yns)| is striking—it is only negative along the approximate phase boundary between positive and negative β regimes, corresponding to β=0.

Combinations of Notch and Jagged production rates lead to random phenotypic decisions (β=0 limit)

In the previous section, we intentionally left out the case of β=0. Interestingly, one can analytically calculate, as shown in the STAR Methods section, that this marginal case is related with the quantity ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|<0. The latter is ratio phenotypic responses when intrinsic over extrinsic variables variations.

Calculating explicitly this aforementioned ratio for different combinations of (D0,J0), we observe that it is negative along the boundary of the two regimes and positive elsewhere (see Figure 7B). This can be interpreted that intrinsic dynamics drive the phenotypic decisions everywhere except the boundary β=0. At the boundary, the cell does not anymore control the decisions/fates and phenotype changes due to microenvironmental fluctuations.

In Figure 8A, we show the different steady-state patterns emerging from different β regimes. We see lateral induction for β >0 cases, resulting in a homogeneous hybrid phase, while, for β <0, we see patches of salt-and-pepper patterned phase surrounded by a sea of cells of the sender phenotype. For β≈ 0, we do not see any discernible pattern. To quantify the deviations of the patterns from a fully random configuration, we calculated the corresponding structure factor, i.e., the power spectrum of the spatial Fourier transform. In turn, we have compared the Area Under the Curve (AUC) between the corresponding pattern of (J0,D0) pair and a purely random pattern (see STAR Methods section). The resulting Figure 8C shows an impressive coincidence with the β sign graph in Figure 8B.Figure 8 Tissue-level pattern formation

The different steady-state configurations of a square lattice of cells correspond to different β values and their population distributions. The color in the heatmaps in the six peripheral plots of (A) represents the level of min-max normalized NICD values ((Xfinal−100)/(1400−100)). The population distributions of NICD values are also captured by kernel density estimates; their modality index (see Equation 22) can be mapped to their β values. (B) is the same as Figure 7A(i), which has a color bar running from −1 to +1. Specifically, the yellow region is β >0 (corresponding to a homogeneous phase of hybrid cells), and the blue region is β <0 (salt-and-pepper islands in a sea of cells of the sender phenotype). For β∼ 0, we see random patterns. This is further underlined by its correspondence with (C), a metric that purely captures the spatial patterning of each population. The metric is the normalized area under the curve of the power spectra of each pattern. The details can be found in Figure S3.

Discussion

In this article, we proposed a mathematical model of lattice-based phenotypic decision-making. We assumed that each cell in the lattice will change its phenotype according to the entropy of its microenvironment. This model is consistent with the recently formulated LEUP theory for cell decision-making in multi-cellular systems.

We obtained a relation for the sensing intensity parameter that could be interpreted to shed light on its possible origin: in the low noise approximation limit, it (the sensing intensity parameter) is the ratio of logarithmic factors of (1) the gradient of sensing dynamics with respect to microenvironmental variables and (2) the variance of the local microenvironmental distribution. Subsequently, we used the example of Notch-Delta-Jagged signaling to estimate the sensitivity parameters. Using this example, we established the conditions for which it could have three possible cases (i.e., β=0, β<0, and β>0). In the agent-based model of phenotypic decisions, we characterized the pattern formations that appeared due to microenvironmental phenotypic interactions and then compared our model with classical mechanistic Notch-Delta-Jagged signaling. To study this, we use two neighborhoods, i.e., (1) the nearest neighborhood (r=1) and (2) the extended neighborhood (r>1) interactions. For the nearest neighborhood interactions, the distribution of phenotypes changes from bimodal to trimodal distributions; as sensing intensity increases, the frequency of hybrid phenotypes increases. We also found chessboard-like patterns that resemble those found in the mechanistic Notch-Delta-Jagged model. Interestingly, for the extended neighborhood (r>1) interactions, we found unimodal distributions as we escalate the sensing intensity values. Additionally, we used the rdfs and average local variance to quantify spatial patterns in the phenotypes and explored the “fluid”-like structures developed in the phenotypic space. Finally, we have quantitatively established that the sign of the β parameter can encode the patterns produced by NJD simulations (quantified via a function of the corresponding structure factor).

The Notch pathway is ubiquitous in nature and is responsible for various biological processes, viz., embryonic development, wound healing, angiogenesis, and cancer progression.6 The intracellular interplay between Notch-Delta-Jagged signaling and transcription factors- and microRNAs-based regulatory network for EMT4 during cancer metastasis might result in other different interesting spatiotemporal patterns at the tissue level. This Notch-EMT-combined mechanistic circuit would have a different microenvironment. However, the Notch-related mechanism is not the sole mechanism involved in EMT regulation, since it can be induced also via Wnt/beta pathway (WNT) and transforming growth factor β pathways.29 To take into account this regulation ambiguity, we decided to study a generalized mechanism that results in the same phenotypic decision. LEUP offers a generalized theory of cell decision, where EMT could be studied without the exact knowledge of the underlying regulatory mechanisms. This paper is the first step toward this goal since proliferation and motility are not included.

Another interesting aspect of our study regards the finding of the “no decision” boundary. In particular, we found out that there is a family of Delta and Jagged production rates where cell decisions are influenced by microenvironmental fluctuations. This constitutes a testable hypothesis since one could design experiments that modulate accordingly the corresponding production rates and observe the corresponding cell fate decisions. If this is true, it can provide ideas for treating pathologies that depend on impaired Notch-Delta signaling. For example, beyond its crucial role in cancer progression, dysregulation of this pathway has also been implicated in various diseases like pulmonary hypertension and Duchenne muscular dystrophy (for a review, see Siebel et al.30). Mutations affecting different components of the pathway have been associated with a spectrum of disorders such as cerebral autosomal dominant arteriopathy with subcortical infarcts and leukoencephalopathy, Alagille syndrome, and spondylocostal dysostosis, among numerous others.30

On a technical note, we have assumed that the distribution of the microenvironment follows the Gaussian distribution. However, one can use a different microenvironmental distribution to study patterns emerging from phenotypic decisions. If a strong coupling exists between the microenvironment and the cell, one may use other definitions of entropy like non-ergodic entropy (e.g., Tsallis entropy,31 Renyi entropy,32 etc.), which will further help us to investigate phenotypic decision-making in a deeper level. Choosing such entropic functionals mandates further theoretical investigations. Finally, another interesting model extension is considering an off-lattice spatial domain to study the effect of LEUP-driven phenotypic decision-making.

An intriguing extension of our research, extending beyond the NDJ pathway, involves mechanisms associated with toggle-switch genes.33,34 These genes are pervasive in developmental processes and also play significant roles in cancer. Toggle-switch mechanisms are characterized by their ability to switch between ON and OFF states, exhibiting bistable behavior similar to that observed in NDJ. Specifically, we calculate our parameter β for various toggle-switch pathways and analyze its relationship to the corresponding phenotypic distribution and pattern formation.

Here, we have not considered the phenomena of cell proliferation or migration. It will be interesting to see how the rate of proliferation can influence the pattern formation regimes. This is particularly interesting for the case of the go-or-grow phenomenon (EMT/Mesenchymal-Epithelial Transition [MET] which have proliferation associated with the epithelial state), which is found in highly invasive cancers such as gliomas.35,36,37,38 Moreover, one can study phenotypic decision-making in different topologies (from networks to lattices and spatial geometries). We have used periodic boundary conditions in the simulations. Further work should explore the impacts of other different boundary conditions (e.g., reflective and absorbing boundary conditions, etc.) on phenotypic decision-making on the spatial scale.

In summary, we have shown how microenvironmental uncertainty can drive the phenotypic decision-making on lattice, which further helps us to study phenotypic decision-making when the precise knowledge about the internal states of the cells is unknown and/or partially known.

Limitations of the study

There exist some limitations of the study, which are as follows. (1) Although the simplification of the entropy-driven principle LEUP is powerful for understanding biological phenomena, it might not fully disclose the multifaceted nature of cellular interactions and signaling pathways. So, one needs biological evidence and knowledge to enhance the theory for further predictions. (2) The exploration of hybrid and trimodal phenotypes is needed in the NDJ model, such that one can encapsulate the complex patterns through LEUP formalism. (3) Lack of experimental pieces of evidence can be an obstacle for this kind of data-driven/information-oriented formalism, which can be tricky in some scenarios.

Resource availability

Lead contact

For additional resources or inquiries please kindly contact Prof. Dr. Haralampos Hatzikirou (haralampos.hatzikirou@ku.ac.ae).

Materials availability

This paper did not include newly generated reagents.

Data and code availability

• The data used in this paper are publicly available.

• All codes and figures are available publicly on the GitHub page of A.A.P. at https://github.com/aditi-pujar/Microenvironmental-Entropy-Dynamics-in-Notch-Delta-Jagged-Signalling-Context.git.

• For further details needed to re-examine the data analysis and methodology presented in this paper, please kindly request the authors for additional information.

Acknowledgments

H.H. and A.B. would like to thank 10.13039/501100001663 VolkswagenStiftung for its support of the “Life?” program (96732 ). H.H. has received funding from the 10.13039/501100002347 Bundesministerium für Bildung und Forschung under grant agreement no. 031L0237C (MiEDGE project/ERACoSysMed). Finally, H.H. acknowledges the support of the RIG-2023-051 grant from 10.13039/501100004070 Khalifa University and the UAE-NIH Collaborative Research grant AJF-NIH-25-KU . M.K.J. was supported by the Ramanujan Fellowship (SB/S2/RJN-049/2018) awarded by the Science and Engineering Research Board, Government of India. U.R. acknowledges the support of the C.V. Raman Fellowship, Indian Institute of Science, Bangalore, and DST-INSPIRE Faculty Fellowship (DST/INSPIRE/04/2022/003052 ), Government of India.

Author contributions

Conceptualization, M.K.J. and H.H.; methodology, A.A.P., A.B., P.S.D., D.S., U.R., and H.H.; software, A.A.P., P.S.D., and D.S.; formal analysis A.A.P., A.B., U.R., P.S.D., and D.S.; investigation, A.A.P., A.B., P.S.D., D.S., M.K.J., and H.H.; resources, M.K.J. and H.H.; writing – original draft preparation, all the authors contributed equally; writing – review and editing, all the authors contributed equally; visualization, A.A.P., D.S., and P.S.D.; supervision, M.K.J. and H.H.; project administration, M.K.J. and H.H.; funding acquisition, M.K.J. and H.H. All authors have read and agreed to the published version of the manuscript.

Declaration of interests

The authors declared no competing interests.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Software and algorithms	
	
Notch-Delta-Jagged signaling mechanistic simulations	This paper	https://github.com/aditi-pujar/Microenvironmental-Entropy-Dynamics-in-Notch-Delta-Jagged-Signalling-Context.git	
LEUP theory driven simulations	This paper	https://github.com/aditi-pujar/Microenvironmental-Entropy-Dynamics-in-Notch-Delta-Jagged-Signalling-Context.git	
Analytical derivations of LEUP theory	This paper	N/A	

Method details

Entropic forces dictate phenotypic dynamics: The LEUP theory

Bayesian cell decision-making

We postulate that cellular decision-making is based on two types of information (see Figure 1) the (1) internal variables (i.e. representing genes, RNA molecules, translational proteins, metabolites, receptors, phenotypes, etc.) that contain the local information sensed by a cell and (2) the external variables (i.e. ligands, chemicals, nutrients, cellular density or stress fields, etc.). In a cell population, for the n-th cell at a time t, the external factors are defined as Ynt∈Rd×N , where N is the number of microenvironmental factors, and internal variables of n-th cell at time t are denoted by the vector Xnt∈Rd. The corresponding dimension is the same because the external variables for a cell can be linear combinations of the same internal variables of the neighbouring cell, e.g. RNA expression. We assume the cell as a Bayesian decision maker who tries to estimate the posterior distribution P(Xnt∣Ynt) at time t:(Equation 10) P(Xnt∣Ynt)=P(Ynt∣Xnt)P(Xnt)P(Ynt)

Here, the local microenvironmental knowledge gained by the cell is represented by the likelihood function P(Ynt∣Xnt). This likelihood function is modified by a multiplicative term of the ratio of the probability distribution of internal variables, P(Xnt), the prior and external variables P(Ynt).

For a cell to construct an accurate likelihood function P(Ynt|Xnt) entails a significant energetic investment by the cell, since its specification requires the production of sufficient membrane receptors. Therefore, defining an informative prior P(Xnt) provides an energetically savvy strategy for cells. A biological example can be the response of the cell in a nutrient-enriched environment. Cells sense their local environment through different biochemical processes, like polymerizing pseudopodia, trans-locating receptor molecules or modifying its cytoskeleton according to mechanical signals.39,40 If cells build an informative prior, then the price for recruitment can be reduced over time t. This precisely the premise of the Least Microenvironmental Uncertainty Principle (LEUP)20,21 which has been linked to the Bayesian Learning in a time independent manner.23

The Least microEnvironmental Uncertainty Principle (LEUP)

In order to calculate the steady state probability distribution function (pdf) of the internal states of the cell, we developed a maximum entropy principle for cellular internal state distributions which will maximize the internal variable entropy constrained by the mutual information I¯(Yns:Xns) between environment and cellular internal variables. The evaluation of the internal states at the steady state in terms of a variational principle can then be written as(Equation 11) δδP(Xns){S(Xns)+β[∫P(Xns)I(Yns:Xns)dXns−I¯(Yns:Xns)]−λ[∫P(Xns)dXns−1]}=0,

The functional derivative with respect to the internal variables is denoted by δ/δP(Xns). The constraints associated with the Lagrange multipliers in Equation 11, i.e., β and λ correspond to the target value of the local mutual information I¯(Yns:Xns) and the normalization constant of the internal state probability distribution. Extra information regarding internal or external variables from biological experiments can be included in terms of statistical observables as additional constraints. Solving Equation 11, the probability distribution of internal states results in a Gibbs-type distribution:(Equation 12) P(Xns)=eβI(Yns:Xns)ZI=e−βS(Yns∣Xns=x)Zs,

where the normalization constant is ZI=∫e−βI(Yns:Xns)dXns and Zs is defined as ∫e−βS(Yns∣Xns)dXns . In SI 3, we present the details of the calculation. The parameter β encodes the sensitivity of the cell to the microenvironment. From Equation 12, we can see that the quantity S(Yns∣Xns) is akin to the potential of the system. In the next section, we shall use this analogy to further write the phenotypic Langevin equation for an arbitrary cell.

Phenotypic Langevin equation

During cellular decision-making, a cell changes its phenotype (internal state) according to the states of its neighbouring cells (microenvironmental state). Through the Langevin equation, one can describe the phenotypic dynamics of the cell, when interacting with its neighbourhood. We assume that the phenotype of the nth cell at time t∈R+ is a continuous variable, i.e., Xnt∈Rd. The external variable, Ynt∈Rd×N is the set of cell phenotypes within a sensing/interaction radius r∈R, i.e. Ynt=(Xn,1t,Xn,2t,…,Xn,Nt), where N is the number of neighbors, as shown in Figure 1B.

Moreover, due to the absence of phenotypic inertia, one can consider an overdamped Langevin equation for the phenotypic time evolution. We consider that the rate of change of the phenotype is associated with the gradient of microenvironmental entropy and Gaussian noise (present during the decision-making). So, the noise has a unit variance and zero mean. The properties of the noise are ⟨ξnX(t)⟩=0 and ⟨ξnX(t1)ξmX(t2)⟩=2Dδ(t1−t2)δmn where t1 and t2 are the two distinct time points. The diffusion constant in the space of phenotypes is denoted by D. The Delta function is defined by δ(t1−t2) and the Kronecker delta by δmn. Recalling the microenvironmental entropy as being analogous to a potential, the Phenotypic Langevin equation for the n-th cell in a population of cells, with phenotype Xnt therefore reads(Equation 13) dXntdt=−β∇XS(Ynt|Xnt)+ξnX(t)

For a small noise approximation, one can consider the above Equation 13 to then take the form of the following ordinary differential equation:(Equation 14) dXntdt=−β∇XS(Ynt|Xnt).

Let us assume that each cell senses/samples the corresponding microenvironmental phenotypes from a Gaussian distribution. The noise term in the Langevin’s equation (i.e., Equation 13) is due to the decision-making stochasticity or to the microenvironmental sensing. From here and onwards, we will assume cells as perfect sensors and we will use the zero noise limit case as represented by Equation 14. Assuming that the phenotypes are all scalars (the number of dimensions of the internal variables, d=1), one can write the phenotype of neighboring cells as the external variables for our n-th cell, i.e. Yn,it=Xn,it,i=1,2,…,N for the N-neighboring cells. Assuming that the microenvironmental entropy term S(Ynt|Xnt)=12log[2πeσYnt|Xnt2], where σYnt|Xnt2 is the sample variance of the cell neighborhood, the entropic force simplifies to:(Equation 15) ∂S(Ynt|Xnt)∂Xnt=1η(XntN−1N2∑jNXjt)

where,(Equation 16) η=Xnt2N+1N∑jXnt2−2XntN2∑jXjt−1N2(∑jXjt)2

In this mathematical model, the key parameters are (i) the sensing intensity parameter β and (ii) the microenvironmental radius of interaction r. The latter defines the size of the sensing neighborhood, i.e., all the cells within a distance r from a cell will contribute to the RHS of Equation 15 and influence the cell’s decision-making. Variations of the two parameters induce a wealth of collective phenotypic dynamics, as shown in the results section.

Mathematical model of Notch-Delta-Jagged signaling

The components of this pathway are the trans-membrane receptor protein, Notch, ligands Jagged and Delta and an intra-cellular factor, NICD (Notch Intra-Cellular Domain). For a given cell, Notch (N) binds to the ligands Delta (D) and Jagged (J) of its neighboring cell. This trans-interaction (with other cells) results in the cleavage of the NICD. NICD then goes to the nucleus and indirectly modulates the production of all three proteins. It transcriptionally activates Jagged and Notch and represses Delta. This asymmetric regulation results in different types of patterns in Delta-dominated versus Jagged-dominated systems. N0,D0,andJ0 represent the basal production rates of the proteins.

The functions Hs+/Hs− represent the transcriptional activation/inhibition (They are shifted Hill functions5 of the form Hs(X)=1+λ[X]n1+[X]n) of the corresponding species due to NICD. The parameter kC corresponds to the cis-inhibition rate. Here, the Notch binds to the Delta or Jagged of the same cell and forms a complex and is hence removed from the system. The parameter kT represents the trans-activation rate which corresponds to the interaction with the neighbouring cell’s Notch, Delta, or Jagged. In the last terms, γ represents the basal degradation rate of the species. We assume first-order kinetics for this.(Equation 17) dNdt=N0Hs+(I,λI,N)−kCN(D+J)−kTN(Dext+Jext)−γNdDdt=D0Hs−(I,λI,D)−kCND−kTDNext−γDdJdt=J0Hs+(I,λI,J)−kCNJ−kTJnNext−γJdIdt=kTN(Dext+Jext)−γII

where Next, Dext, Jext refers to the number of Notch, Delta, and Jagged proteins on the parts of the surfaces (that are in contact with the cell) of the neighbouring cells. λI,i (i=N, D, J) denotes the maximum fold changes in the expression of the respective quantities due to NICD.

Quantification of microenvironmental entropy/variance in Notch-Delta-Jagged system

Let us formalise the notion of local microenvironmental entropy by defining an order parameter as ⟨σ2⟩L. We define ⟨σ2⟩L to be the average local variance for a population of cells. For our n-th cell with Xnt=Xnt (the one dimensional phenotype variable) and Ynt=(Yn1t,Yn2t,…,YnNt)=(Xn1t,Xn2t,…,XnNt) (external variables are the phenotypes of its N neighbours), the mean external phenotype it senses is Xnt¯. And their variance is,(Equation 18) σn2=Var(Ynt)=∑iN(Xnit−Xnt¯)2N,⟨σ2⟩L=1|L|(∑nεLσn2).

where L refers to the set of cell indices on the lattice and the corresponding cardinality |L| is the number of lattice points. Thus, ⟨σ2⟩L is the local variance, averaged over the entire lattice.

This is a measure of local entropy. This arose from postulating that in LEUP, the conditional distributions are normally distributed, implying that their entropy is given by a function of the variance.(Equation 19) S(Ynt|Xnt=Xnt)=12log(cσ2(Ynt|Xnt=Xnt))=12log(cσn2)

(Equation 20) S¯(Ynt|Xnt)=12log(c⟨σ2⟩L)

(Equation 21) ⇒⟨σ2⟩L=(eS¯(Ynt|Xnt))2/c

where c=2πed/2. Thus, our average local variance actually is a measure of the constrained conditional entropy, S¯(Ynt|Xnt) at a given point in time.

Quantification of distribution multi-modality

The difference between lateral inhibition and lateral induction shows up in the internal state distribution of cell populations; those with Delta-dominant signalling (lateral inhibition) are bimodal whereas those with Jagged-dominant signalling (lateral induction) are unimodal.

Sarle’s Bimodality Index41 is a statistical index between 0 and 1 that distinguishes between unimodal, uniformly random and multimodal distributions:1. SBC>0.55¯ for multimodal distributions.

2. SBC=0.55¯ for uniformly random distributions.

3. SBC<0.55¯ for unimodal distributions.

Therefore, we concretely define our Modality Index as(Equation 22) ModalityIndex=0.55¯−SBC.

Simulation details

Using the Notch-Delta-Jagged mechanistic model and LEUP driven phenotypic Langevin equation, we simulate the ODEs that govern both Notch Delta and LEUP evolution on cells on fixed square lattices (of size 24×24). The time steps taken were dt=0.01s and the simulations were run for 50,000 to 300,000 iterations.

For the LEUP simulations, each cell was assigned a random variable Xnt∈[−1,1] that captured the phenotype of the cell. The LEUP phenotypes then evolved in accordance with our phenotypic Langevin equation, (13) for different parameter regimes (varying interaction radius r and sensitivity β). In the Notch Delta simulations, each cell was initialised with levels of Notch, Delta, Jagged, and NICD for different basal production rates of Delta and Jagged (i.e. D0 and J0) and evolved over time. While all the proteins were important markers, we chose the NICD level to be the variable that captured the overall state of the cell best. In other words, the Notch Delta Jagged phenotype of a cell was its NICD concentration. Please note that in all simulations we use periodic boundary conditions.

Formulating LEUP in terms of constrained mutual information

Recalling the variational Equation 11 and defining,I¯(Xns,Yns)=S(Yns)−S¯(Yns|Xns)⇒I¯(Xns,Yns)=S(Yns)−S¯(Yns|Xns)

The constraint on conditional entropy can be equivalently formulated as a constraint on mutual information. Thus our variational equation becomes,(Equation 23) δδP(Xns){S(Xns)+β[∫P(Xns)(S(Yns)−S(Yns∣Xns))dXns+(S(Yns)−S¯(Yns∣Xns))]}=0

Which simplifies to,δδP(Xns)[S(Xns)+β(∫dXnsP(Xns=x)I(Yns|Xns=x)−I¯(Yns|Xns))]=0

Which implies P(Xns=x)=eβI(Yns|Xns=x)ZI. To conclude,(Equation 24) ∴P(Xns=Xnt)=eβI(Yns,Xns=Xnt)ZI=e−βS(Yns|Xns=Xnt)ZS

where,(Equation 25) ZI=∫e+βI(Yns,Xns=x)dx=eβS(Yns)∫e−βS(Yns|Xns=x)dx=eβS(Yns)ZS

Thermodynamic identity

Previously we have derived:P(Xnt)=eβI(Yns,Xns=x)ZI⇒lnP(Xnt)=βI(Yns,Xns=x)−lnZI

Now,S(Xns)=−∫P(Xns=x)lnP(Xns=x)dx=lnZI∫P(Xns=x)dx−β∫P(Xns=x)I(Yns,Xns=x)dx

Therefore, we can derive the corresponding thermodynamic identity for the phenotypes s:(Equation 26) S(Xns)=lnZI−βI(Yns,Xns)

Time evolution of ⟨σ2⟩L

In Figure S1, we provide detailed simulations for the temporal evolution of ⟨σ2⟩L for different parametric regimes.

Sensitivity parameter β calculation

By definition, we have,Z′=∫|ΣYns∣Xns|−β/2dXns=(2πe)βd/2∫((2πe)d|ΣYns∣Xns|)−β/2dXns=(2πe)βd/2∫e−βS(Yns|Xns)dXns

From the identity S(X,Y)=S(X)+S(Y|X)≤S(X)+S(Y), we know,(Equation 27) S(Y|X)≤S(Y)

We now have to deal with positive and negative β separately.

The β>0 case

Continuing from Equation 27,βS(Yns|Xns)≤βS(Yns)⇒∫e−βS(Yns|Xns)dXns≥∫e−βS(Yns)dXnsZ′≥(2πe)βd/2e−βS(Yns)(∫dXns)

Setting M=∫dXns, where M is some characteristic concentration of internal variables, we ultimately arrive at the inequality,(Equation 28) Z′≥Zbound′

(Equation 29) whereZbound′=|ΣYns|−β/2M

therefore, for β>0, we obtain a lower bound on its partition function. From Equation 8, we have,β2|ΣYns∣Xns|=ln|∇YG(Xns,Yns)|+S(Yns)−d/2−lnZ′⇒β2|ΣYns∣Xns|≤ln|∇YG(Xns,Yns)|+S(Yns)−d/2−lnZbound′⇒β(|ΣYns∣Xns|−|ΣYns|2)≤ln|∇YG(Xns,Yns)|+S(Yns)−c⇒β(−I(Xns,Yns))≤ln|∇YG(Xns,Yns)|+S(Yns)−c∴β≥−(ln|∇YG(Xns,Yns)|+S(Yns)−cI(Xns,Yns))

We use the property that I(Xns,Yns)≥0 in the last step. Thus, the estimate for β we obtain is actually a lower bound for β>0.

The β<0 case

We can work out the above calculation in an exactly analogous way in order to arrive at,Z′≤Zbound′β≤−(ln|∇YG(Xns,Yns)|+S(Yns)−cI(Xns,Yns))

i.e. Zbound′ is an upper bound for our partition function, and the estimate we derive for β is actually an upper bound for β<0.

Thus, by substituting in Zbound′ for Z′, we arrive at lower bounds for β>0 and upper bounds for β<0. In other words, we arrive at the lower bounds for |β| in both regimes. In the main text, we are going to use limiting value for all values of β.

The β formula for NDJ system

In this section, we shall consider a system of equations representing Notch-Delta-Jagged signaling5,27 with the presence of Gaussian noise. We denote the Notch expression of the n-th cell at time t by Nnt, its Delta expression by Dnt, its Jagged expression by Jnt and its NICD expression by Int. We have seen how the system evolves in the ODEs outlined in the Mathematical Framework. We define two new variables here, Lnt=Dnt+Jnt and Lext=Dext+Jext, the total internal and external ligand concentrations. We obtain how Lnt evolves over time by adding the ODEs for Dnt and Jnt. Thus, our new set of variables - Notch Nnt, Ligand (Lnt) and NICD (Int) levels evolve according to the following ODEs -(Equation 30) dNntdt=N0Hs+(I,λI,N)−kCNntLnt−kTNntLext−γNnt+ξndLntdt=D0Hs−(I,λI,D)+J0Hs+(I,λI,J)−kCNntLnt−kTLntNext−γLnt+ξndIntdt=kTNntLext−γIInt

In the adiabatic approximation, the internal state of the cell changes at a much faster rate than the average of its neighbours. Hence, the external variables are taken as time-independent in Equation 30. The constants N0 and D0 are the production rates of Notch and Delta proteins, which are measured in the order of molecules/h. I0 is the threshold value of the NICD complex Hill function in a number of molecules.5,27

As the timescale of NICD molecule formation happens on a faster time scale one can easily assume that the rate of change of NICD molecules reaches a steady state (faster than Notch, Delta, and Jagged molecules). In the quasi-steady state, the amount of NICD complex molecule can be written as,(Equation 31) Int=kTγINntLext=αNntLext

The parameter α=kTγI. This approximation holds because NICD degrades rapidly as compared to N, D and J.42 Now, we can substitute the value of Int in Equation 30 to get a simplified equation for Notch-Delta-Jagged signaling as,(Equation 32) dNntdt=N0Hs+(αNntLext,λI,N)−kCNntLnt−kTNntLext−γNnt+ξndLntdt=D0Hs−(αNntLext,λI,D)+J0Hs+(αNntLext,λI,J)−kCNntLnt−kTLntNext−γLnt+ξn

We expand the first term (the shifted Hill function) in the Taylor series of Equation 32 around the dissociation constant, I0. We specifically did the Taylor series expansion around the dissociation constant to understand the switching behaviour of internal states and the correlation between internal and environmental variables. Considering the first two terms in the Taylor series expansion, we can further write Equation 32 approximately as,(Equation 33) dNntdt=(λI,N+1)N02+hN0(λI,N−1)4I0(αNntLext−I0)−kCNntLnt−kTNntLext−γNnt+ξndLntdt=(λI,D+1)D02+hD0(λI,D−1)4I0(αNntLext−I0)+(λI,J+1)J02+hJ0(λI,J−1)4I0(αNntLext−I0)−kCNntLnt−kTLntNext−γLnt+ξn

or in the matrix form, the above system of equations can be written as,(Equation 34) ddt(NntLnt)=((λI,N+1)N02+hN0(λI,N−1)4I0(αNntLext−I0)…−kCNntLnt−kTNntLext(λI,D+1)D02+hD0(λI,D−1)4I0(αNntLext−I0)…+(λI,J+1)J02…+hJ0(λI,J−1)4I0(αNntLext−I0)−kCNntLnt−kTLntNext)−γ(NntLnt)+(ξnξn)

Absorbing γ into the characteristic time scale of the system, we can now define the internal states as,(Equation 35) Xnt=(NntLnt),

and external states as,(Equation 36) Ynt=(NextLext)

Using this definition of internal and external variables, we can write the system of Equation 34 as,(Equation 37) ddtXnt=G(Ynt,Xnt)−Xnt+ξ¯

where,(Equation 38) G=((λI,N+1)N02+hN0(λI,N−1)4I0(αNntLext−I0)…−kCNntLnt−kTNntLext(λI,D+1)D02+hD0(λI,D−1)4I0(αNntLext−I0)…+(λI,J+1)J02…+hJ0(λI,J−1)4I0(αNntLext−I0)−kCNntLnt−kTLntNext)=(G1G2)ξ¯=(ξnξn)

To find the steady state condition, we take dNntdt=dLntdt=0. Now, we can use the tools previously mentioned in Equation 2 as,(Equation 39) Xns=G(Yns,Xns)+ξ¯

Finding the drift term ∇YG(Xns,Yns)

In order to calculate β, the dynamics-related term we need to calculate is ∇YG(Xns,Yns), (which takes on the meaning of drift.)∇YG(Xns,Yns)=|∂G1∂Next∂G1∂Lext∂G2∂Next∂G2∂Lext|=|0Nns(hN04I0α(λI,N−1)−kT)−kTLnsNnshD04I0α(λI,D−1)+NnshJ04I0α(λI,J−1)|=kT2(N0h(λI,N−1)4I0γI−1)|NnsLns|

This finally simplifies to,(Equation 40) ∇YG(Xns,Yns)=|kTκ|NnsLns||

where, κ=kT(1−hNN0(λN−1)4I0). (This is a convenient way to define these parameters because of a later calculation involving ∇XG(Xns,Yns)). The values of the parameters we used have been derived from5,27 and have been listed in the S.I.

Numerical calculation of β from NDJ simulations

We now proceed to perform the above calculation for a given point in the Notch Delta Jagged phase space. Recalling the (fully expanded) β formula,β=−(ln|∇YG(Xns,Yns)|+S(Yns)−(lnM+d/2)12(ln|ΣYns|−ln|ΣYns∣Xns|))

We obtained the formula for the drift term, ∇YG(Xns,Yns) in the previous subsection. There are a few more statistical quantities we need to calculate. We obtain them thus - In our Notch Delta Jagged systems (note that d=2),Xns=(NL)andYns=(NextLext)

Quadratic term correction in β calculation

There is initially an additional term of the form (Yns−F(Yns))TΣY−1(Yns−F(Yns)) in the numerator of β that we assumed went to zero. We calculated it as a correction term - F(Yns) was estimated to be the average of a cells’ neighbours’ Yns’s and the term hence calculated. We obtained values of the order 10−17 (indeed, comparable to machine precision). We can therefore conclude that the assumption to neglect this term is well grounded.

Calculating ΣYns

This statistical term appears in the denominator as well as the numerator (as S(Yns)=12ln((2πe)2|ΣYns|)). By the above definition of the Yns variable, we have,(Equation 41) ΣYns=(Cov(Next,Next)Cov(Next,Lext)Cov(Lext,Next)Cov(Lext,Lext))

where for each entry of the matrix, the covariances of two variables are calculated over the whole lattice. |ΣYns| is simply the determinant of the above matrix.

Calculating ΣYns∣Xns

This is a conditional covariance matrix and its computation is slightly more involved. For multivariate normal distributions - which we have assumed throughout our treatment of all these variables, the conditional covariance matrix is given by the following identity -,(Equation 42) ΣYns∣Xns=ΣYns−ΣYnsXns(ΣXns)−1ΣXnsYns

For our Notch Delta Jagged system,(Equation 43) ΣXns=(Cov(N,N)Cov(N,L)Cov(L,N)Cov(L,L))

(Equation 44) ΣXnsYns=(Cov(N,Next)Cov(N,Lext)Cov(L,Next)Cov(L,Lext))

(Equation 45) ΣYnsXns=(Cov(Next,N)Cov(Next,L)Cov(Lext,N)Cov(Lext,L))

Thus, we can find |ΣYns∣Xns| by taking the determinant of the above conditional covariance matrix.

Calculating c

This is a trivial calculation - we know d=2 here and the characteristic levels of Notch and the ligands can be taken to be 5000 units/cell. Thus,c=lnM+d/2=ln(500012+12)+2/2=≃7.071×103

Having thus obtained the methods for calculating each of the terms in our formula for β, we are now in a position to calculate β over sweeps of the Notch Delta Jagged landscape.

Phase space of other terms in β formula

In Figure S2, we present how the individual terms that contribute to the analytical β calculation are quantified in the (D0,J0) phase space.

Calculation details for β=0 condition

In the following, we detail the calculations for deriving the condition for β=0:ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|<0.

Approximation of the internal variable variance

Consider a dynamical system where a cell with internal variables Xnt in a microenvironment with external variables Ynt is subject to the following ODE,(Equation 46) Xnt.=G(Xnt,Ynt)+ζ(t)

Taking Taylor expansion around the stationary point (Xns,Yns), where the system is at equilibrium,(Equation 47)

Taking squares and integrating on both sides,(Equation 48) ∫(Xnt−Xns)2P(Xnt,Ynt)dXntdYnt=|∇YG(Xns,Yns)∇XG(Xns,Yns)|2∫(Ynt−Xns)2P(Xnt,Ynt)dXntdYntσX2=|∇YG(Xns,Yns)∇XG(Xns,Yns)|2σY2.

Upper bound for entropy, S(X) from dynamics

Now, we know the mean and variance of our cell’s internal variables, Xns from Equations 47 and 48. We can therefore assume P(Xns) to be the maximally entropic probability distribution for the given constraints, which is the normal distribution. i.e.(Equation 49) P(Xns)∼N(Xns−∇YG(Xns,Yns)∇XG(Xns,Yns)(Ynt−Yns),|∇YG(Xns,Yns)∇XG(Xns,Yns)|2σY2)⇒S(Xns)≤12ln((2πe)dσXns2)S(Xns)≤12ln((2πe)d|∇YG(Xns,Yns)∇XG(Xns,Yns)|2σYns2)

where, d is the dimension of the variable Xnt and σYns2 is the variance of Yns. Here, we make the implicit assumption that Ynt is also normally distributed. Thus, its variance, σYns2, corresponds to that of a multivariate normal distribution. i.e. σYns2=|ΣYns|, the determinant of its covariance matrix.

Phase space boundary calculation

We know that the boundary between the positive and negative β regimes would be given by β=0. Thus, to theoretically ground the observation from Figure 7, we are required to show that,β=0⇒ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|<0

At the phase boundary, we are given β=0. Recalling Equation 26,(Equation 50) ⇒S(Xns)=lnZI≤S(Yns)+ln|∇YG(Xns,Yns)∇XG(Xns,Yns)|⇒ln|∇YG(Xns,Yns)∇XG(Xns,Yns)|≥lnZI−S(Yns)∴ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|≤S(Yns)−lnZI

Forβ=0

lnZI=∫eβI(Yns,Xns=x)dXns=∫dXns=M

where M is some characteristic number of molecules ∼103. Thus,ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|≤S(Yns)−M

We can see from the phase portrait of log|ΣYns| that,(Equation 51) log|ΣYns|∼(−8,25)

(Equation 52) ⇒ln|ΣYns|∼(−82.303,252.303)

(Equation 53) ⇒S(Yns)=12ln((2πe)2|ΣYns|)∼(1.10101,8.26558)

(Equation 54) ∴ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|≦S(Yns)−M≪0

Therefore,(Equation 55) β=0⇒ln|∇XG(Xns,Yns)∇YG(Xns,Yns)|<0

We could arrive at the phase boundary result theoretically/empirically as well.

Spatial patterning quantification via structure factor

To quantify the correspondence between population distributions (captured by the modality index) and spatial patterning, we follow strategy based on the spatial Fourier transform (FT). We first obtain the two dimensional FTs of the spatial pattern. These are plotted in Figure S3B.IFFT(k)=∑x∑yINDJ(x,y)e−2πi(kxx+kyy)

We then obtain the amplitude spectrum - by averaging all the FT values for a given wavenumber k=(kx,ky). (All the lattice points at a given radial distance, equidistant from the origin of the Fourier transform lattice correspond to the same magnitude of |k|). The amplitude spectrum when squared gives us the power spectrum, shown in Figure S3A.pκ(k)=∑x∑yIFFT(k)2δ(kx2+ky2−k)

For a given pattern, we integrate the area under the curve of its power spectrum, AUCNDJ and “normalise” it by subtracting the area under the curve of a random reference spatial pattern AUCrandom.AUCNDJ=∑kpκ(k)

We see in Figure S3A that while Delta dominated chessboard patterns have a larger area under the curve (corresponding to β<0) and Jagged dominated patterns have smaller areas (corresponding to β>0), along the phase boundary (where β≈ 0, the AUCNDJ is almost equal to AUCrandom. Thus, the phase diagram of this metric, AUCNDJ−AUCrandom which captures the spatial patterning, has an almost one-to-one correspondence with the phase diagram of both, modality index and calculated β values (Figure 8).

Parameter values for Notch Delta Jagged and LEUP simulations

Notch Delta Jagged simulations

These are taken with reference to,41. Basal transcription rates:

N0=1400D0∈[500,2000]J0∈[200,1300]

It is also mentioned that the largest number of Notch proteins as well as Jagged/Delta molecules in a cell are ∼5000. Thus,M=|Xns|=50002≈7×103

2. Rate constants:

γ=0.1γI=0.5

kT=5×10−5kC=5×10−4

3. Shifted Hill function parameters:

ThresholdvalueofNICD,I0=200Co−operativity,nN=2.0nD=2.0nJ=5.0Foldchange,λN=2.0λJ=2.0λD=0.0

Note that a shifted Hill function is of the form,Hs(I)=1+λ(II0)n1+(II0)n

4. Lattice parameters:

SizeofLattice=25×25dt=0.01sT=300sNumberofiterations=300/0.01=30000

LEUP simulations,SizeofLattice=25×25dt=0.01sT=50sNumberofiterations=50/0.01=5000

Quantification and statistical analysis

No statistical tests were performed in the study.

Supplemental information

Document S1. Figures S1–S3 and Table S1

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110569.
==== Refs
References

1 Zhou B. Lin W. Long Y. Yang Y. Zhang H. Wu K. Chu Q. Notch signaling pathway: Architecture, disease, and therapeutics Signal Transduct. Targeted Ther. 7 2022 95 10.1038/s41392-022-00934-y
2 Tasca A. Helmstädter M. Brislinger M.M. Haas M. Mitchell B. Walentek P. Notch signaling induces either apoptosis or cell fate change in multiciliated cells during mucociliary tissue remodeling Dev. Cell 56 2021 525 539.e6 10.1016/j.devcel.2020.12.005 33400913
3 VanDussen K.L. Carulli A.J. Keeley T.M. Patel S.R. Puthoff B.J. Magness S.T. Tran I.T. Maillard I. Siebel C. Kolterud Å. Notch signaling modulates proliferation and differentiation of intestinal crypt base columnar stem cells Development 139 2012 488 497 10.1242/dev.070763 22190634
4 Boareto M. Jolly M.K. Goldman A. Pietilä M. Mani S.A. Sengupta S. Ben-Jacob E. Levine H. Onuchic Jose’N. Notch-Jagged signalling can give rise to clusters of cells exhibiting a hybrid epithelial/mesenchymal phenotype J. R. Soc. Interface 13 2016 20151106 10.1098/rsif.2015.1106 27170649
5 Jolly M.K. Boareto M. Lu M. Onuchic J.’N. Clementi C. Ben-Jacob E. Operating principles of Notch–Delta–Jagged module of cell–cell communication New J. Phys. 17 2015 055021 10.1088/1367-2630/17/5/055021
6 Bocci F. Nelson Onuchic J. Jolly M.K. Understanding the Principles of Pattern Formation Driven by Notch Signaling by Integrating Experiments and Theoretical Models Front. Physiol. 11 2020 929 10.3389/fphys.2020.00929 32848867
7 Reppas A.I. Lolas G. Deutsch A. Hatzikirou H. The Extrinsic Noise Effect on Lateral Inhibition Differentiation Waves ACM Trans. Model. Comput. Simul. 26 2016 1 18 10.1145/2832908
8 Yasugi T. Sato M. Mathematical modeling of Notch dynamics in Drosophila neural development Fly 16 2022 24 36 10.1080/19336934.2021.1953363 34609265
9 Koon Y.L. Zhang S. Rahmat M.B. Koh C.G. Chiam K.-H. Enhanced Delta-Notch lateral inhibition model incorporating intracellular Notch heterogeneity and tension-dependent rate of Delta-Notch binding that reproduces sprouting angiogenesis patterns Sci. Rep. 8 2018 9519 10.1038/s41598-018-27645-1 29934586
10 Xu X. Seymour P.A. Sneppen K. Trusina A. Egeskov-Madsen A.la R. Jørgensen M.C. Jensen M.H. Serup P. Jag1-Notch cis-interaction determines cell fate segregation in pancreatic development Nat. Commun. 14 2023 348 10.1038/s41467-023-35963-w 36681690
11 Sprinzak D. Lakhanpal A. LeBon L. Santat L.A. Fontes M.E. Anderson G.A. Garcia-Ojalvo J. Elowitz M.B. Cis-interactions between Notch and Delta generate mutually exclusive signalling states Nature 465 2010 86 90 10.1038/nature08959 20418862
12 Koizumi Y. Iwasa Y. Hirashima T. Mathematical study of the role of Delta/Notch lateral inhibition during primary branching of Drosophila trachea development Biophys. J. 103 2012 2549 2559 10.1016/j.bpj.2012.11.005 23260057
13 Kang T.-Y. Bocci F. Nie Q. Onuchic J.N. Levchenko A. Spatial–temporal order–disorder transition in angiogenic Notch signaling controls cell fate specification Elife 12 2024 RP89262 10.7554/eLife.89262.3
14 Cohen R. Taiber S. Loza O. Kasirer S. Woland S. Sprinzak D. Precise alternating cellular pattern in the inner ear by coordinated hopping intercalations and delaminations Sci. Adv. 9 2023 eadd2157 10.7554/eLife.89262.3 36812313
15 Balázsi G. van Oudenaarden A. Collins J.J. Cellular Decision Making and Biological Noise: From Microbes to Mammals Cell 144 2011 910 925 10.1016/j.cell.2011.01.030 21414483
16 Gigante G. Giuliani A. Mattia M. A novel network approach to multiscale biological regulation Cell Syst. 14 2023 177 179 10.1016/j.cels.2023.02.004 36924765
17 Rørth P. Fellow travellers: emergent properties of collective cell migration EMBO Rep. 13 2012 984 991 10.1038/embor.2012.149 23059978
18 Smart M. Zilman A. Emergent properties of collective gene-expression patterns in multicellular systems Cell Rep. Phys. Sci. 4 2023 101247 10.1016/j.xcrp.2023.101247
19 Hatzikirou H. Statistical mechanics of cell decision-making: the cell migration force distribution J. Mech. Behav. Mater. 27 2018 1 7 10.1515/jmbm-2018-0001
20 Barua A. Nava-Sedeño J.M. Meyer-Hermann M. Hatzikirou H. A least microenvironmental uncertainty principle (LEUP) as a generative model of collective cell migration mechanisms Sci. Rep. 10 2020 22371 10.1038/s41598-020-79119-y
21 Barua A. Syga S. Mascheroni P. Kavallaris N. Meyer-Hermann M. Deutsch A. Hatzikirou H. Entropy-driven cell decision-making predicts ‘fluid-to-solid’ transition in multicellular systems New J. Phys. 22 2020 123034 10.1088/1367-2630/abcb2e
22 Barua A. Beygi A. Hatzikirou H. Close to Optimal Cell Sensing Ensures the Robustness of Tissue Differentiation Process: The Avian Photoreceptor Mosaic Case Entropy 23 2021 867 4300 10.3390/e23070867 34356408
23 Barua A. Hatzikirou H. Cell decision-making through the lens of Bayesian learning Preprint at arXiv 2023 10.48550/ARXIV.2301.06941
24 Friston K. Is the free-energy principle neurocentric? Nat. Rev. Neurosci. 11 2010 605 10.1038/nrn2787-c2
25 Knill D.C. Pouget A. The Bayesian brain: the role of uncertainty in neural coding and computation Trends Neurosci. 27 2004 712 719 10.1016/j.tins.2004.10.007 15541511
26 Bialek W. Biophysics: Searching for Principles 2012 Princeton University Press
27 Boareto M. Jolly M.K. Lu M. Onuchic J.N. Clementi C. Ben-Jacob E. Jagged–Delta asymmetry in Notch signaling can give rise to a Sender/Receiver hybrid phenotype Proc. Natl. Acad. Sci. USA 112 2015 E402 E409 10.1073/pnas.1416287112 25605936
28 de Back W. Zhou J.X. Brusch L. On the role of lateral stabilization during early patterning in the pancreas J. R. Soc. Interface 10 2013 20120766 10.1098/rsif.2012.0766
29 Dongre A. Weinberg R.A. New insights into the mechanisms of epithelial–mesenchymal transition and implications for cancer Nat. Rev. Mol. Cell Biol. 20 2018 69 84 10.1038/s41580-018-0080-4
30 Siebel C. Lendahl U. Notch signaling in development, tissue homeostasis, and disease Physiol. Rev. 97 2017 1235 1294 10.1152/physrev.00005.2017 28794168
31 Tsallis C. Possible generalization of Boltzmann-Gibbs statistics J. Stat. Phys. 52 1988 479 487 10.1007/BF01016429
32 Buckland W.R. Neyman J. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability. Volume 1: Theory of Statisticts Statistician 13 1963 159 10.2307/2986976
33 Giuliani A. Bui T.T. Helmy M. Selvarajoo K. Identifying toggle genes from transcriptome-wide scatter: A new perspective for biological regulation Genomics 114 2022 215 228 10.1016/j.ygeno.2021.11.027 34843905
34 Alon U. Network motifs: theory and experimental approaches Nat. Rev. Genet. 8 2007 450 461 10.1038/nrg2102 17510665
35 Alfonso J.C.L. Köhn-Luque A. Stylianopoulos T. Feuerhake F. Deutsch A. Hatzikirou H. Why one-size-fits-all vaso-modulatory interventions fail to control glioma invasion: in silico insights Sci. Rep. 6 2016 37283 10.1038/srep37283
36 Alfonso J.C.L. Talkenberger K. Seifert M. Klink B. Hawkins-Daarud A. Swanson K.R. Hatzikirou H. Deutsch A. The biology and mathematical modelling of glioma invasion: a review J. R. Soc. Interface 14 2017 20170490 10.1098/rsif.2017.0490
37 Hatzikirou H. Basanta D. Simon M. Schaller K. Deutsch A. ‘Go or grow’: the key to the emergence of invasion in tumour progression? Math. Med. Biol. 29 2012 49 65 10.1093/imammb/dqq011 20610469
38 Mascheroni P. López Alfonso J.C. Kalli M. Stylianopoulos T. Meyer-Hermann M. Hatzikirou H. On the impact of chemo-mechanically induced phenotypic transitions in gliomas Cancers 11 2019 716 10.3390/cancers11050716 31137643
39 Sourjik V. Berg H.C. Receptor sensitivity in bacterial chemotaxis Proc. Natl. Acad. Sci. USA 99 2002 123 127 10.1073/pnas.011589998 11742065
40 Ueda M. Shibata T. Stochastic Signal Processing and Transduction in Chemotactic Response of Eukaryotic Cells Biophys. J. 93 2007 11 20 10.1529/biophysj.106.100263 17416630
41 Thomas R.K. Bimodality revisited J. Mod. Appl. Stat. Methods 6 2007 3 10.22237/jmasm/1177992120
42 Wang M.M. Notch signaling and Notch signaling modifiers Int. J. Biochem. Cell Biol. 43 2011 1550 1562 10.1016/j.biocel.2011.08.005 21854867
