Energy landscape decomposition for cell differentiation with proliferation effect

Abstract Complex interactions between genes determine the development and differentiation of cells. We establish a landscape theory for cell differentiation with proliferation effect, in which the developmental process is modeled as a stochastic dynamical system with a birth-death term. We find that two different energy landscapes, denoted U and V, collectively contribute to the establishment of non-equilibrium steady differentiation. The potential U is known as the energy landscape leading to the steady distribution, whose metastable states stand for cell types, while V indicates the differentiation direction from pluripotent to differentiated cells. This interpretation of cell differentiation is different from the previous landscape theory without the proliferation effect. We propose feasible numerical methods and a mean-field approximation for constructing landscapes U and V. Successful applications to typical biological models demonstrate the energy landscape decomposition's validity and reveal biological insights into the considered processes.


I. THEORY OF THE ENERGY LANDSCAPE DECOMPOSITION
We begin with the cell population equation where x ∈ R n is the gene expression vector, b(x) is the drift term, ϵ is the noise amplitude, R(x) is the birth and death rate (BDR), and p(x, t) is the probability density function (pdf) of cells over the gene expression space. The energy landscape decomposition (ELD) is to find energy landscapes U (x), V (x), and curl part term f (x) which satisfy The decomposition has infinite choices, but a meaningful one is supposed to use U to describe the cell types and use V to describe the differentiation. There are data-based approaches and model-based approaches for constructing energy landscape. In this study, ELD is a model-based approach with assumption that b(x), ϵ and R(x) are known, which means we can solve Eq. S1 or simulate a corresponding stochastic system Eq. S10 to get the final steady pdf, denoted as P U (x). Thus, U is defined by U can be taken as the simplest gradient system whose equilibrium density equals to P U (x). Metastable states on U stand for different cell types. On the other hand, by setting R(x) as zero, we can obtain another steady pdf denoted as P 0 (x). V is defined by V can be taken as the change of landscape U caused by the cell birth and death. Values of V can represent the pluripotency, and the differentiation evolves from the high to the low value of V . The remained part in b(x) is the curl part term defined as f satisfies the property ∇ · (f (x)P 0 (x)) = 0, which corresponds to the divergence free condition of the curl flux J (x) = f (x)P 0 (x) defined in [1,2].

A. Gradient system
A special case of Eq. S1 is b(x) = −∇F (x), which is a gradient system. We can obtain the analytic expression of P 0 (x) as where Z 0 is a normalization constant. By Eq. S3, Eq. S4, Eq. S5, and Eq. S6, we can obtain F (x) = U (x) + V (x), in the sense of ignoring a constant, and f (x) = 0. Further, through the Eq. S3 and the Eq. S1 at t = +∞, we can obtain that the potential U satisfies a forward equation and the potential V satisfies another backward equation Thus, for a gradient system, the decomposition Eq. S2 degenerates to the methods used in the population balance analysis (PBA) [3] and the landscape of differentiation dynamics (LDD) [4]. The difference is that PBA and LDD are data-based, and estimate U and V by solving Eq. S7 and Eq. S8 from real datasets. The EDL in this paper is model-based and tries to estimate U and V by Eq. S3 and Eq. S4. As a remark, for a non-gradient system, Eq. S7 always holds, while Eq. S8 doesn't. Through we know Eq. S8 holds if and only if f ⊥ ∇V . Estimating V by solving Eq. S8 in PBA and LDD can not be generalized to the non-gradient system.

B. No birth and death system
When the system has no birth and death, i.e. R(x) ≡ 0, we can obtain P 0 (x) ≡ P U (x) and V (x) ≡ 0. Thus the decomposition Eq. S2 reduces to the landscape proposed by Wang et al. [1,2]. But Wang's landscape is constructed for each set of possible parameters, such as different strength of gene-gene interaction or cell volume. The change of the landscape U during the differentiation is controlled by the manually setting parameters. However, the parameters in our ELD can be fixed. If the essential proliferation and death of cells are taken into account, i.e. R(x) ̸ ≡ 0, we show that an induced landscape V can characterize the direction of differentiation.

C. Weighted stochastic equation form
There exists a weighted stochastic dynamics corresponding to the population equation Eq. S1 as where X t (y, ω) is a trajectory ω with the initial point y, and ρ t (y, ω) is the weight for this trajectory. Each trajectory represents a cell in the differentiation process. All the trajectories have the same weight as t = 0. If the weight ρ of ω increases, it is supposed that the cell has proliferation, while if the weight ρ of ω decreases, death is supposed to occur. The population pdf is computed as the weighted expectation of all trajectories, i.e. p(x, t) = E {ρ t (y, ω)δ(x − X t (y, ω))}, where δ(·) is the Dirac delta function.

Proposition:
The population pdf of stochastic system Eq. S10 satisfies Eq. S1. Proof: By Ito formula, we have which is the same as Eq. S1. The initial p(x, 0) is set as the distribution of y. □

II. NUMERICAL ALGORITHM FOR THE LOW DIMENSIONAL SYSTEM
The key point to obtain the decomposition is to estimate the two steady pdfs P U (x) and P 0 (x). For the low-dimensional system, we can use numerical algorithms, such as finite difference or finite element method, to solve Eq. S1. In the following, we will show the details of streamline diffusion method [5], which is applied in our examples.
The weak form of Eq. S1 in streamline diffusion method is where ⟨·, ·⟩ is the integration over the considered spatial space Ω, v + δv b is the test function, v ∈ H 1 0 (Ω), v b = b · ∇v, and δ is a small constant. In our 2-dim examples, we use the bilinear function over rectangle grid in the spatial axis. On each grid element, we assume b(x), ∇ · b(x), R(x) are constants with the value at the left-bottom point. If we choose the implicit 1-order discretization of the temporal axis and denote the basis over spatial grid as {φ i (x)}, the Eq. S12 is discretized as where τ is the time step, and p t i is the value of p(x, t) at the grid point x i . Eq. S13 is a linear equation, which can be denoted as , and N is the number of grid nodes in Ω. For each square finite element K with spatial step size h, the bilinear bases are With Dirichlet boundary condition, the global A and B can be assembled by the local matrices and where After sufficient iterations of Eq. S14 with normalized p at each time step, we can obtain an estimation of the steady distribution p(x) on the grid.
There are two remarks: 1. The purpose of the streamline diffusion method is to reduce the oscillation caused by the numerical algorithm, especially around the shock wave domain. 2. δ is usually chosen asch, where h is the spatial grid size,c is a small positive number if ϵ < h, or 0 if ϵ ⩾ h.

III. MEAN-FIELD APPROXIMATION IN THE HIGH DIMENSIONAL CASE
It is not possible to solve Eq. S1 directly in the high-dimensional case. Instead, a mean-field approximation (MFA) can give a computable estimation of the steady pdfs P U and P 0 as a Gaussian mixture form where K is the number of components obtained by the counts of stable states in the deterministic dynamicsẊ t = b(X t ), is the kth Gaussian component with mean µ (k) and covariance ϵΣ (k) , and ρ (k) is the mixture weight with K k=1 ρ (k) = 1.

A. MFA for the weighted dynamical system
In this subsection, we show the MFA for Eq. S10, which is also applied in our numerical example. There are two parts in the MFA Eq. S19. Gaussian component Eq. S20 is obtained by short time asymptotics, while mixture weight ρ (k) is obtained by long time asymptotics.

Short time asymptotics
In O(1) time scale, the particle/cell stays in one attractive basin without transition. We can have the approximation where (x − µ t ) (k) is a k-dimensional tensor whose j = (j 1 , j 2 , . . . , j k )-th element is k s=1 (x js − µ js ), and the operator A : B = j A j B j is the summation of corresponding element-wise product of two tensors. We can find that X t is independent of R(x) in Eq. S10. Then the evolution equations for µ t and Σ t are and where δ ij is the Kronecker delta function with value 1 when i = j and 0 when i ̸ = j. Through different initial points, we can obtain all the components (µ (k) , Σ (k) ) by finding the stable points of the determined system where L 1 = (L i 1 ) n×1 is the evolving operator for µ t , and L 2 = (L ij 2 ) n×n is the evolving operator for Σ t .

Long time asymptotics
To determine the mixture weights {ρ (k) }, we need to consider the transitions between attractive basins in the long time scale log(t) ≳ O(1/ϵ). In this regime, X t is upscaled to a continuous-time Markov chain, which transits between different basins. Denote the Taylor expansion of R(x) as then the average BDR for state k is Let Q be the upscaled generator matrix, in which Q ij is the transition rate from state i to state j when i ̸ = j, and Q ii = − j̸ =i Q ij . We can have the asymptotics for mixture weights as U ) contains the mixture weights for P U (x), and ρ 0 = (ρ (k) 0 ) contains the mixture weights for P 0 (x).
The accurate Q is difficult to be obtained, so we use the approximation where Ω k denotes the k-th stable state, and X 0 (ω) is sampled under a uniform distribution on a finite domain. The approximation Eq. S29 can be taken as the percentage of trajectories falling into the stable state k (size of the basin).
To estimate the weight ρ i . Set q k = −Q kk as the exit rate of state k, then h i obeys the exponential distribution with rate parameter q k .
Then, by ergodic theorem we can obtain where M is the number of T i k in [0, T ]. To ensure the non-explosion of the probability, we must assume that lim t→∞ where h is the random variable standing for the holding time, and q k ⩾ R k to ensure the nonexplosion of pdf. Further, we assume the exit time τ k for each component in F (x) is in proportion to ρ (k) 0 . By The constant C can be obtained by optimizing which is induced from the steady constraint R(x)P U (x) dx = 0. ρ U , the steady pdfs P 0 (x) and P U (x) can be estimated by Eq. S19. Further, landscapes U (x) and V (x) are obtained by Eq. S3 and Eq. S4, respectively.

B. MFA for the population-level equation
There is MFA directly applied to the population-level equation Eq. S1, too. The short time asymptotics for the Gaussian components is and For the Gaussian case, the forth order moments can be calculated by The estimations for weights are and where X t (ω) follows Eq. S10, Ω 0 k is the kth state for P 0 , and Ω U k is the kth state for P U . However, µ t and Σ t depends on R(x) now. When the system is with or without R(x), the approximated global pdf Eq. S19 may have quite different number and position of the components (i.e. P 0 (x) and P U (x) can have disparate components). The severe consequence is that the energy landscape V (x) could have a very bad performance to describe the differentiation after so many procedures with approximation. That's the reason why we propose MFA to Eq. S10 for multi-state systems in the applications.
But, for the system with only one state/component, Eq. S34 and Eq. S35 could differ P 0 and P U , while the MFA to the weighted dynamics Eq. S25 fails.

C. Analytic solution of the mean-field approximation for two examples
In the MFA, we expand the dynamics b(x) and R(x) to the order O(ϵ), instead of O(1) which is used in reference [2]. We will use two analytic examples in this subsection to show its necessity, in order to get an accurate approximation to the landscape V (x). An intuitive illustration in the gradient system is given by Eq. S8, which says for an example, we obtain V (x) ∼ O(ϵ) and the approximation to the dynamics of P 0 and P U needs to be at least O(ϵ).

1-dim gradient system
In the first example, we use a one dimensional gradient system with true potentials The task is to decompose the system with known (S40) We can find that only µ ∞ = 0, Σ ∞ = 1 is the stable solution, and the approximated steady pdf is where Z U is a normalization constant. The estimated energy landscape U by Eq. S3 is U (x) = 1 2 x 2 + ϵ log(Z U ), which equals to U (x) in the sense of ignoring a constant. The mean-field approximation to ∂p(x,t) We can get the approximated steady pdf P 0 (x) = 1 Thus, in this example, we show that the mean-field approximation can get the correct decomposition of the system for a one-dimensional Ornstein-Uhlenback process with source and sink/birth and death. As a remark, if the approximations in Eq. S34 and Eq. S35 are not in the order of O(ϵ) but only O(1), the estimation of P U in Eq. S40 will be incorrect and we may obtain spurious U (x), V (x).

2-dim non-gradient system
In the second example, we use a two dimensional non-gradient system with true terms U (x) = The task is to decompose the system with known By solving Eq. S34 and Eq. S35 at t = +∞, we can obtain a stable point µ ∞ = [0, 0] T and Σ ∞ = I when the non-zero R(x) exists, where I is the 2-dim unit matrix. Thus, the mean-field approximation of the steady pdf is P U (x) = 1 2πϵ exp − 1 2ϵ (x 2 1 + x 2 2 ) . The energy landscape U is estimated by Eq. S3 as U (x) = 1 2 (x 2 1 +x 2 2 )+ϵ log(2πϵ), which equals to U (x) in the sense of ignoring a constant.
When setting R(x) ≡ 0, by solving Eq. S34 and Eq. S35 at t = +∞, we can obtain a stable point µ ∞ = [0, 0] T and Σ ∞ = 1 1+ϵ I. The mean-field approximation of the steady pdf is P 0 ( In this example, we show that the mean-field approximation can get the correct decomposition of the system for a two-dimensional non-gradient process with source and sink/birth and death. It can be also found that the expansions to the order O(ϵ) in Eq. S34 and Eq. S35 are necessary for the correct estimation of U , V and f .

A. Drift-diffusion process
The dynamics of the 2-dim drift-diffusion process is which is a gradient system with potential F (x) = , and white noise with amplitude √ 2ϵ can be added to the system. ϵ = 0.01. BDR is set as

B. Two-gene fate decision
The dynamics for the two-gene fate decision system is where the two genes x 1 and x 2 have self-activation and inter-inhibition, white noise with amplitude √ 2ϵ can be added to the system, and parameters are set as α 1 = α 2 = 0.3, β 1 = β 2 = 0.5, n = 4, k 1 = k 2 = 1, and S = 0.5. ϵ = 0.01. BDR is set as where the rate amplitude r can vary from 0 to 30.

C. T-cell differentiation process
The dynamics for the T-cell differentiation [6] is modeled as The Hill function x n x n +K n represents for activation and K n x n +K n stands for inhibition. Notch signal α N , which regulates x 1 , x 3 and x 4 , is fixed as zero in our simulation. ϵ = 0.001. BDR is set as . We use the mean-field approximation Eq. S25 to calculate the ELD. In order to simplify the computing complexity, we only consider diagonal elements of Σ t , and Eq. S25 can be reduced to dρ t (ω) = R(X t (ω))ρ t (ω) dt, ∂ jj R(µ (k) ) · Σ jj (k) , when X t (ω) in state k. (S46)

D. 10-dimensional system with the birth and death term
We use 10-dim dynamics to show that the MFA is applicable to the ELD in high-dimensional space with where X(t) ∈ R 10 , the potential function is noise amplitude is ϵ = 10 −4 , and W (t) is the standard Brownian motion. There are three metastable states around (1, 1, . . . , 1) T , (−1, −1, . . . , −1) T , and (0, 0, . . . , 0) T , respectively. We use a BDR term When we use the MFA Eq. S25 to calculate the ELD, the results of landscapes U and V are shown in Fig. S2. U can indicate the metastable states, while V can imply the pluripotency or differentiation direction successfully.