Cellular gradient flow structure connects single-cell-level rules and population-level dynamics

In multicellular systems, the single-cell behaviors should be coordinated consistently with the overall population dynamics and functions. However, the interrelation between single-cell rules and the population-level goal is still elusive. In this work, we reveal that these two levels are naturally connected via a gradient flow structure of the heterogeneous cellular population and that biologically prevalent single-cell rules such as unidirectional type-switching and hierarchical order in types emerge from this structure. We also demonstrate the gradient flow structure in a standard model of the T-cell immune response. This theoretical framework works as a basis for understanding multicellular dynamics and functions.


INTRODUCTION
Multicellular systems are organized dynamically and robustly to shape various populational patterns required to achieve biological functions [1][2][3].Because the population dynamics is realized by single-cell-level processes, i.e., cellular proliferation, death, migration, and differentiation, behaviors of individual cells should be coordinated consistently with the overall population dynamics, which may rule the single-cell processes.
For example, phenotypic switching and differentiation of a cell in a population are often unidirectional.Moreover, the multiple cell types are hierarchically ordered, and their kinetic properties, e.g., type-switching and proliferation rates, also seem to be coupled.The existence of hierarchy, or equivalently, acyclic cell-type lineage structures and kinetic coupling are prevalent in multicellular systems from immunity to development [4][5][6][7].These single-cell rules emerging in a population may be related to the functions and the coordination of the population.
The notion of the epigenetic landscape has been used pervasively to describe the directional dynamics among hierarchical cellular types, where individual cells are likened to balls rolling down the landscape [8][9][10].However, the interrelation of such single-cell-level dynamics and landscape with the overall population dynamics and their function has been rarely investigated except in a few prescient works pointing out its importance [11][12][13].
In this work, we show that the interrelation between single-cell and population levels can emerge from a gradient flow structure of a population.We model the desirable population distribution for achieving a biological function by the landscape of a utility function.Then, we derive the population dynamics that maximizes the utility given biological costs of single-cell processes, which results in a gradient flow of the utility function.We demonstrate as an example that the standard model of Tcell population dynamics in the acute immune response [14] can be understood as the gradient flow.From the populational gradient flow structure, the single-cell-level landscape emerges, from which the unidirectional type switching, hierarchical cell-type, and kinetic couplings are generally derived.Moreover, the single-cell landscape is related to the population-level utility landscape as its functional variation.Our result can work as a theoretical basis to bridge single-cell-level rules and behaviors with population-level dynamics and functions of multicellular systems.

GRADIENT FLOW OF CELLULAR POPULATION DYNAMICS
We firstly introduce the governing equations of the heterogeneous cellular population dynamics.We consider a large heterogeneous population of cells with different types.The state of population at time t is characterized by n t = {n t (x)} x∈X , where n t (x) ≥ 0 is the population size of the cells with type x ∈ X.Here X is a set of all possible types.The state of the population changes over time by the following cellular actions: growth (proliferation minus death), type switching, and immigration (recruitment) of new cells from outside of the population.Their rates are assumed to be type-dependent such that g t (x) is the growth rate of type x, v t (x, y) ≥ 0 is the type-switching rate from type x to type y, and m t (x) ≥ 0 is the immigration rate of type x.In this work, we focus on the case where X is discrete, but it can be easily extended to the continuum case.Then, the population dynamics of the cells can be described by the following equations where we approximate n t (x) as a continuous variable, which is valid if the population size is large enough.Next, we connect the population dynamics with biological functions.To this end, we introduce a utility arXiv:2205.13143v1[q-bio.PE] 26 May 2022 function U t (n).The utility function U t (n) abstractly represents how good a given population distribution n is under the situation at time t.The time dependence of U t (n) is essential for modeling various biological situations.For example, if a pathogen invades our body, a particular population distribution n of immune cells would work more effectively than another n , which is represented as U t (n) > U t (n ).The utility may change after the eviction of the pathogen, which is captured by the time dependence of U t (n).Another example of the utility function is for developmental processes.A specific pattern of differentiated cells would be required at time t, which may change as the development progresses.Therefore, the dynamics that can induce the population distribution n t (x) into the one with a higher utility more promptly would be more functional than other dynamics.
However, the rates of growth, immigration, and type switching cannot be arbitrarily high due to the biological cost of those processes and physical constraints.In order to account for it, we introduce the cost function C n (g, m, v).The cost function abstractly characterizes the instantaneous biological cost of taking cellular actions at given rates (g, m, v), which is nonnegative and dependent on the current population n.In this work, we consider the cost function to have the following form where w g , w m , and w v are positive weights whose values depend on the single-cell level mechanisms of actions.This form of the cost function is derived from three assumptions: 1) costs for different cellular actions are independent, i.e., the total cost is just a sum of them, 2) growth costs and type-switching costs are proportional to the current cell number, and 3) for each cellular action, the cost is a strictly convex smooth function of the rate.The second assumption is reasonable because growth and type-switching costs are incurred for each cell in the current population.In contrast, since immigration is usually independent of the current population, the immigration cost is not proportional to the current cell number.The third assumption is crucial to prohibit the optimal action rates from being arbitrarily high.Such unrealistic behavior can happen if the cost grows slowly as the rates increase, and the utility can cancel it out.While we focus here on the quadratic cost function, the simplest convex function, our theory can be extended to more general convex functions.
We consider the population dynamics of Eq. ( 1) where the rates (g t , m t , v t ) are determined to maximize the value of the utility function under the cost.
where Diff nt U t (g, m, v) is the time derivative of U t through the time evolution of n t given rates (g, m, v): We could obtain the explicit form of the unique optimum (g † t , m † t , v † t ) of the above optimization problem as follows (see Appendix for derivation) where [a] + is the positive part of a ∈ R, and ∇ is the discrete gradient operator, i.e., for any φ : X → R, ∇φ(x, y) := φ(y) − φ(x).Note that these optimum rates are scale-invariant: they are invariant under the rescaling of the utility function U t and the cost function C with the same factor.Let us consider the dynamics with the optimal rates (g t , m t , v t ) = (g † t , m † t , v † t ).Under this dynamics, the value of the utility function U t (n) evolve as Since the instantaneous cost C nt is nonnegative, the value of the utility function always increases when ∂Ut ∂t vanishes.Indeed, when U t does not depend on t, the dynamics is a generalized gradient flow of U , where C is called as a dissipation function [15][16][17].

T-CELL IMMUNE RESPONSE MODEL
We demonstrate that a model of T-cell population dynamics in the acute immune response [14] is a gradient flow in our sense.Here we introduce a slightly modified version of the model proposed in [14].The model assumes three types of T cells, naive (N ), activated effector (A), and memory (M ), and the numbers of these cells n t (N ), n t (A), n t (M ) are described by the following ordinary differential equations where g, m, v are constant growth, immigration and typeswitching rates, and I t represents the temporal change in the environmental situation such that I t = 0 when the immune cells should contract to recover to the normal state, and I t = 1 when they should expand to eliminate pathogens (Fig. 1 (a)).We assume on/off transition: A typical time evolution is depicted in Fig. 1 (b).In the expansion phase (τ 0 ≤ t < τ 1 ), the number of activated T cells rapidly increases, whereas, in the contraction phase (τ 1 ≤ t), the number of memory T cells increases instead.
There are several versions of the T-cell immune response model, and their mathematical properties were investigated in [18][19][20].The model we introduced in Eqs. ( 6)-( 7) is simple yet includes immunologically realistic factors.Notably, a simplified version of the model introduced here was shown to reproduce experimental data [14,21].In Appendix, we list some of the other Tcell immune response models and discuss their gradient flow structures.
The model Eqs.( 6)-( 7) is a gradient flow in our framework on three-type space X = {N, A, M } with the following utility function U t and the cost function C. The utility function is a time-dependent linear function with coefficients u N := g N ρ N w 0 , u (It) A := (g A0 + g A1 I t )w 0 , and u M := g M ρ M w 0 , indicating that each cell type has different importance depending on the environmental situation.The cost function (Eq.( 2)) is specified with the weights and all the other weights are +∞.Here, w 0 , ρ N , and ρ M are arbitrary positive constants.By taking account of the scale invariance of the optimal rates, w 0 is the scaling factor for the utility function and the cost function.We define it as the same as the growth weight for activated T cells.ρ N and ρ M are the relative growth weights for naive and memory T cells.
To demonstrate that the model is actually a gradient flow, we numerically calculated U t (n t ), the value of utility function along the time evolution of n t (Fig. 1 (c)).The result shows that U t (n t ) is monotonically increasing except at the change point t = τ 1 , where ∂Ut ∂t in Eq. ( 5) becomes nonzero.Thus, the gradient flow structure completely explains this behavior.
Finally, we discuss that the gradient flow structure is qualitatively consistent with biologically plausible parameters.The structure imposes the positivity of the weights (Eqs.( 9)) of the cost function.Rearranging the terms, we obtain the following constraints on the growth  3) are allowed in the gradient flow, graphs ( 2) and ( 4) contain a cycle and never appear in the gradient flow.For graph ( 2) and ( 4), one of the edges are incompatible with the ordering by δU t (n t ) δn rates: If these inequalities do not hold, the increase in utility is no longer guaranteed.The existence of expansion and contraction in the immune response naturally requires that the growth rate for the activated T cells in the contraction phase be negative (g A0 < 0) and the growth rate for the activated T cells in the expansion phase be positive (g A0 + g A1 > 0).If, in addition, the growth rate g N for the naive T cells is positive, there exists the weight parameters ρ N and ρ M satisfying the constraints.According to [22], naive T cells in humans have a relatively high proliferation rate to maintain the size of the naive Tcell population, implying the growth rate g N for the naive T cells to be positive.Thus, the gradient flow structure is consistent with biologically plausible parameters.

EMERGING SINGLE-CELL RULES
While our gradient flow is derived from the optimization at the population level, i.e., what kind of populational change is better than others, it also determines the behaviors of each cell, i.e., what a cell of type x should do or should not do.We will show two such rules derived from the gradient flow structure.In both rules, the functional derivative δUt(nt) δn (x) of the utility function plays a vital role.One can interpret this function as a utility function at the single-cell level because it defines which type is preferable to the other types.It is in contrast to the original utility function U t (n), which defines a metric only for the population of cells.

Unidirectional phenotypic transition
In multicellular systems, unidirectional phenotypic switchings are commonly observed.For example, the Tcell immune response model does not have bidirectional or cyclic type switchings.We show that such unidirectionality is tightly linked to the gradient flow structure.
To this end, we consider a type-switching graph G t .The nodes of the graph G t are types X, and an edge from type x to y exists if and only if the type-switching rate is not zero, i.e., v † t (x, y) > 0. One can show that this graph G t for any gradient flow is always acyclic (Fig. 2).In the simplest case, the types x and y cannot have bidirectional type switching because v † t (x, y) and v † t (y, x) cannot be simultaneously positive from Eq. (4c).Depending on the sign of ∇ δUt (nt)  δn (x, y), either or both of v † t (x, y) and v † t (y, x) is zero.In other words, cells of type x switch to type y if and only if y is better than x, which means that If we place all the types X vertically in the order of δUt(nt) δn as in Fig. 2, every type-switching edge points upward.Thus, more generally, three or more types cannot have any cyclic type switching.If a cycle exists, at least one of the edges points downward, which contradicts the ordering by δUt(nt)  δn .One can view δUt(nt) δn as a kind of epigenetic landscape of single cells [8].We note that it depends on time t and the current population n t , which could be interpreted as variations of the epigenetic landscape due to timedependent external signals and cell-cell interactions.

Coupling
Growth, immigration, and phenotypic switching in multicellular systems are not independent.For example, in the T-cell immune response model, the growth and type-switching rates change simultaneously when the environmental condition changes.We show that the gradient flow structure implies cooperative relationships among growth, immigration, and type-switching rates.
Consider the simplest setting where all the weights are finite and constant, w g (x) = w m (x) = w v (x, y) = 1 ∀x, y ∈ X.From the fact that the growth rates, immigration rates, and type-switching rates (Eqs.( 4)) have the same term δUt(nt) δn , we find some relations among them.When there is an immigration flux to type x, the growth rate of type x must be positive and vice versa: When the type-switching rate from type x to y is positive, the growth rates of the source type x must be lower than the growth rate of the destination type y One can intuitively understand these effects as cooperation among growth, immigration, and type-switching to achieve the same goal: maximizing the utility function.
The T-cell immune response model has this cooperative coupling as long as the inequalities in Eq. ( 10) are satisfied.The coupling between growth and typeswitching was also predicted in Furusawa and Kaneko's model: the growth rate of stem-type cells is lower than differentiated-type cells [12].Moreover, this kind of coupling has been observed in cell biology: cells undergoing differentiation stop the cell cycle and do not divide [7].Other coupling properties predicted from the gradient flow structure can be used to search for the structure in actual biological systems.

DISCUSSION
Unveiling a potential gradient flow structure of a given population dynamics is reduced to identifying the utility and cost functions.It would be desirable to have a systematic way of identifying these functions.An approach is to infer the utility function from experimental data by using techniques in machine learning and single-cell omics [23,24].Another approach is to derive the cost function from the physical and thermodynamic principles, e.g., by the large deviation theory [17].In either case, our framework will serve as a basis for linking the single-cell and population properties.
Finally, it should be noted that a given population dynamics may not always fall into the class of gradient flow in the strict sense.Some modifications of the T-cell model (Eqs.( 6)-( 7)) can violate the conditions to be a gradient flow.Nevertheless, the gradient-flow-like behaviors can still be preserved if the modification is moderate, and the utility monotonically increases in time (see Appendix for more detail).Thus, our theory can be used to search for such behaviors.Moreover, we can further extend the notion of gradient flow [25] to accommodate oscillatory components, e.g.cell cycle, and others.It expands the applicability of our approach to a wide range of multicellular phenomena and will be pursued.

FIG. 2 .
FIG.2.Examples of type-switching graphs of three types X = {A, B, C}.While acyclic graphs (1) and (3) are allowed in the gradient flow, graphs (2) and (4) contain a cycle and never appear in the gradient flow.For graph (2) and (4), one of the edges are incompatible with the ordering by δU t (n t )