Cell lineage branching as a strategy for proliferative control

Background How tissue and organ sizes are specified is one of the great unsolved mysteries in biology. Experiments and mathematical modeling implicate feedback control of cell lineage progression, but a broad understanding of what lineage feedback accomplishes is lacking. Results By exploring the possible effects of various biologically relevant disturbances on the dynamic and steady state behaviors of stem cell lineages, we find that the simplest and most frequently studied form of lineage feedback - which we term renewal control - suffers from several serious drawbacks. These reflect fundamental performance limits dictated by universal conservation-type laws, and are independent of parameter choice. Here we show that introducing lineage branches can circumvent all such limitations, permitting effective attenuation of a wide range of perturbations. The type of feedback that achieves such performance - which we term fate control - involves promotion of lineage branching at the expense of both renewal and (primary) differentiation. We discuss the evidence that feedback of just this type occurs in vivo, and plays a role in tissue growth control. Conclusions Regulated lineage branching is an effective strategy for dealing with disturbances in stem cell systems. The existence of this strategy provides a dynamics-based justification for feedback control of cell fate in vivo. See commentary article: http://dx.doi.org/10.1186/s12915-015-0123-7. Electronic supplementary material The online version of this article (doi:10.1186/s12915-015-0122-8) contains supplementary material, which is available to authorized users.


Disturbance modeling
All the disturbances and parameter perturbations discussed in the main text are modeled by modifying the right-hand side of (1) to include external disturbance ∆(t)=[∆ 1 (t), ∆ 2 (t)] where ̅ 1 , ̅ 2 are constants; σ 1 (t), σ 2 (t), can be stochastic or deterministic processes; and δ Dirac (t) is the Dirac delta function.
To facilitate the study of the system's response to disturbances, we start by decomposing the disturbance into a time-independent static component ̅ =[ ̅ 1 , ̅ 2 ] and a time-dependent component δ(t)=[δ 1 (t), δ 2 (t)], so that (see Figure B1). The choice of ̅ is not unique and is selected here so that the fluctuations of ∆(t) are centered around ̅ . We are interested in how the population of terminal cells x 2 changes as a function of Δ (t) and how different choices of p r can help attenuate/reject such disturbances -i.e., minimize the perturbations of x 2 for the nominal (undisturbed, Δ(t)=0 ) steady state value ̅ 2 * . If we define ̅ 2 to be the steady state response of the terminal cell population x 2 (t) only to the static component of the disturbance, then we can decompose the disturbance response as follows: i.e., the response of the system to the disturbance is also be decomposed into a static response ̅ 2 =K pr ( ̅ ) and a dynamical response ξ 2 (t)=H pr (t,∆,x(0)) ( Figure B1).
K pr is a (static) nonlinear function of that characterizes the dependence of the steady state value x 2 as a function of a static disturbance ( Figure S1A), x 2 = K pr ( ̅ ), and can be easily computed by solving algebraic equations: ( )  Figure S1).
H pr is a dynamic nonlinear function of the disturbance that characterizes the behavior of the system around the steady state as a function of the disturbances. For a specific choice of parameters and disturbances one

Dynamics near steady state
Under the assumption of small disturbance δ(t), the dynamics of ξ 2 (t) are approximated by , which can be rewritten as: The dynamical systems for both fate and renewal control models can be more instructive in the Laplace (frequency) domain, where it also reveals some non-trivial properties of the system. Indeed, the effect of the disturbances δ 1 , δ 2 on ξ 2 (the perturbations of population of terminal cells from its nominal value) for renewal control is given by Consider the case when the disturbance is a sinusoid δ 1 (t) = sinω 1 t. Then ||ξ 2 || ∞ : = sup t ξ 2 (t), the maximum perturbation of the terminal cells from the nominal value x 2 , is given by |G 1 (jω 1 )| where j = √ , i.e., the (S7) magnitude of G 1 at the frequency ω 1 . Similarly if δ 2 (t) = sinω 2 t, ||ξ 2 || ∞ = |G 2 (jω 2 )|. In general, any disturbance can be expressed as a sum of sinusoidals and the sinusoidal frequencies ω compose what is called the disturbance spectrum. So for a general disturbance, the reduction of the perturbations ||ξ 2 || ∞ is dependent on the magnitude of G 1 and G 2 across the disturbance spectrum. Let ∥G∥ ∞ : = max jω |G(jω)|.
Then ∥G 1 ∥ ∞ and ∥G 2 ∥ ∞ are the maximum perturbation from the nominal value for any frequency of δ 1 and δ 2 respectively.
W 1 and W 2 are independent of p r , and therefore the question of which choice of p r best rejects/attenuates a disturbance δ is equivalent to which choice of p r best reduces the magnitude of S across the disturbance spectrum. We show that for the branched topology, there are constrains that limit how small S can be made.

Renewal control model
In this case the specific transfer functions are given by If L(z) = 0 then z is called the zeros of the system and if S(λ) = ∞, λ is called a pole of the closed loop system. In this case z = v > 0 is a zero of the system, and the system is called non-minimum phase (since z > 0 is a right half-plane (RHP) zero). The steady state error to a small constant perturbation disturbance δ 1 is given by |ξ 2 /δ 1 |=G 1 (0)=h 1 ( ̅ 1 , ̅ 2 )/(2αd ̅ 2 ), i.e., it is inversely proportional to the gain α. Therefore, more aggressive controllers result in smaller errors for disturbances of this type. In general we want G 1 and G 2 to be small across all frequencies ω (not just ω = 0). However the existence of the RHP zero z = v, imposes some hard constraints on how small G 1 and G 2 can be made. First, using maximum modulus theorem can be shown that and therefore, there is a limit to how small the peaks of G 1 and G 2 can be made, independent of the choice of controller.
Furthermore, S must satisfy a special form of Bode's integral formula Since the weights W 1 and W 2 are independent of the controller, then (S8) is a general constraint on G 1 and G 2 (ln|G i (jω)| = ln|W i (jω)| + ln|S(jω)|, i = 1, 2). |S(jω)| < 1 implies an attenuation of noise at frequency ω, while |S(jω)| > 1 implies and amplification of noise at frequency ω. Equation (S8) is a type of conservation law, stating that the net disturbance attenuation and amplification must be balanced. The weight 2v in (S8) is a low pass filter, which makes the disturbance attenuation at low frequencies more costly (i.e., higher amplification at other frequencies). (S8)

Fate control model
In this case the specific transfer functions are given by ( , )

Regulated cell cycle rate for unbranched lineages
Consider renewal control, where there is some additional regulation of the rate at which cells divide, i.e., v is a function of x 2 . The dynamics given by ( ) Let v: = v( ̅ 2 ) and β = − ∂v( ̅ 2 )∂ ̅ 2 . For a given β, the dynamics of (S9) near the steady state are given by

 
The response to these disturbances is given by the transfer functions G 1 = W 1 S and G 2 = W 2 S (for δ 1 and δ 2 respectively), where

Trans-differentiation or delayed differentiation as a "temporary" alternate fate
Here we consider a branched lineage in which the alternative differentiated cell fate x 3 , at some later time, reconnects back to the original differentiated fate, either through trans-differentiation or delayed differentiation so that ultimately there is only one "terminally" differentiated cell type (x2). Might this scheme reap the benefits of fate control without incurring the extra "cost" of producing a potentially unnecessary alternative cell type?
The ODE model for such a scheme would be given by ( ) where u is the rate at which x 3 differentiate (or trans-differentiate) into x 2 (we can add an additional death rate term into the third equation, but it doesn't change the substance of the argument that follows). If we apply the fate control scheme described above, with p d =κp r and negative feedback on p r , we find that as u gets very small, these equations approach the basic fate control model, whereas as u gets very large, the third equation becomes increasingly irrelevant, and the system approaches the renewal control model.
Shown in the left panel is the response for a single sample path realization of the stochastic process. The mean and the standard deviation for 5 different realizations are shown in the right panel. Supporting Tables   Table S1. Parameter values used in the simulations shown in Figures 2 and 4.
The dynamics are given by Eq. (S10). The disturbances enter the system at t=0. For simulations in Figure 5, at t=0 the system is at the steady state x 1 =7.5, x 2 =100, x 3 =7.5, x 4 =75, x 5 =5, x 6 =60 and η (i) (t) , i=1,2,3 are birth-death processes with birth rate 2.5 and death rate 0.025 . For simulations in Figure 6E-6F, the initial population of stem cells x 1 is 10 and there are no other types of cells. The desired final concentrations are 60000, 45000, and 30000 for x 2 , x 4 , and x 6 respectively.  Figure 4) p d (x 2 )=0.5 p r (x 2 ). The disturbances enter the system at t=0 . At t=0 the system is at steady state x 1 (0)=10, x 2 (0)=100. η 1 (t) is a birth-death process with birth rate 2.5 and death rate 0.025. η 2 (t) is a birth-death process with birth rate 2.5 and death rate 0.25. DS1, DS2, DS5 and DS6 are disturbance types described in the Supporting Information. Figure S1