A Minimal Transcriptional Controlling Network of Regulatory T Cell

Subscriber access provided by MONASH UNIVERSITY

Article

A Minimal Transcriptional Controlling Network of Regulatory T Cell Development Chen Liao, and Ting Lu J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp402306g • Publication Date (Web): 03 May 2013 Downloaded from http://pubs.acs.org on May 7, 2013

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A Minimal Transcriptional Controlling Network of Regulatory T Cell Development Chen Liao and Ting Lu∗ Department of Bioengineering and Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA E-mail: [email protected]

Phone: +1 217-333-4627. Fax: +1 217-265-0246

∗ To

whom correspondence should be addressed

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Regulatory T cells (Treg) are a subpopulation of T cells that are central to immune homeostasis and develop under the control of a complex regulatory network consisting of FoxP3 and its partner factors. A central question about this network is how does it enable T cells to robustly specify and stably maintain their states despite intrinsic and environmental fluctuations. Inspired by recent experimental advances, we propose here a minimal transcriptional controlling network and use it to illustrate the robustness and dynamic features of Treg development. Our study shows that the controlling network may exhibit distinct dynamics depending on its parameter regimes and that the maintenance of multistability requires the orchestration of both its positive and negative feedback loops. In addition, system volume contributes monotonically to the increase of the network’s robustness. We further show that the dynamics of our model varies with the alteration of FoxP3-DNA binding affinity, consistent with recent experimental findings. This minimal model thereby offers new insights into the dynamics and robustness of Treg development and may serve as a platform for future exploration towards a more quantitative and systematic understanding of the immune system.

Keywords: Treg signature; functional robustness; network dynamics; effective potential landscape; feedback loops; stochasticity.


Page 2 of 32

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Regulatory T cells (Treg) are a subpopulation of T cells that modulate the immune system, maintain tolerance to self-antigens, and prevent autoimmune disease. 1–3 A fundamental question relating to Treg is how do they develop to implement their vital functionality 4 ? Well established studies have revealed that FoxP3, the characteristic transcription factor expressed by Treg, is central to the establishment and maintenance of immunological self-tolerance and homeostasis , 5–7 suggesting that FoxP3 may serve as an unique window to uncover the development of Treg lineage. FoxP3 was initially considered the ‘master regulator’ of Treg, as supported by studies showing that FoxP3 overexpression enables normal CD4+ T cells to exhibit Treg-like phenotypes while its mutation resulting in the defective development of functional Treg. 8–11 However, recent experiments, which combine chromatin immune precipitation, genome-wide tiling-array analysis and next-generation sequencing, have suggested that FoxP3 only accounts for partial Treg signature. 12–19 Interestingly, one latest study showed that FoxP3 forms large heterogeneous complexes with other partner molecules, including transcriptional factors and chromatin regulators, to achieve promotive or repressive gene regulation, indicating that FoxP3 and its partners constitute a regulatory network possessing both positive and negative feedback loops. 20 In parallel, Fu et al. identified five transcription factors that interact with FoxP3 and act together to form a transcriptional network consisting of multiple and redundant feedback loops. 21 These two studies collectively suggest that Treg development is not determined solely by individual regulatory components but a complicated network of FoxP3 and its interacting partners which self-enhance and function together for developmental control. 22 Such a network picture of the Treg controlling program raises a set of questions relating to network architecture and its associated biological functionality, among which those regarding network robustness is particularly fascinating. As demonstrated by several experimental studies, the Treg lineage is remarkably stable under different physiologic conditions and with various challenges: These CD4+ T cells maintain robustly in either a conventional state (called Tconv) or a regulatory state (called Treg) unless triggered by signals. 23–26 From a dynamical systems perspective, 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

it means that the controlling program of the Treg lineage operates very robustly against intrinsic molecular fluctuations and external environmental perturbations. Why is the regulatory network so robust? What are the architectural characteristics of the network that give rise to this robustness? And what are the dynamic features of the network? Although few computational efforts have been made towards the understanding of T-helper and induced regulatory T cells, 27–30 system-level exploration of the robustness of the controlling program of Treg development has been missing, primarily due to the lack of experimental knowledge until very recently. 20,21 Robustness is the ability of a system to resist change under perturbation or unusual conditions, which is employed as a canonical quantity for characterizing a dynamical system and is also often associated with stochasticity. 31–37 Stochasticity is a ubiquitous feature of biological networks, including those for transcriptional control in gene regulation. 38–40 Stochasticity arises from the random nature of biochemical events engaged in biological networks as well as the inherent bulk variability of cellular environment. 41 As a result, it profoundly influences cellular behavior and development. For instance, noise can perturb a system from its optimal steady state and even possibly drive cells to dramatically different phenotypic states. 42–45 To reveal the consequences of fluctuations on cellular functionality, many theoretical tools, such as those from statistical mechanics and nonlinear dynamics, 46–50 and computational tools, including stochastic simulation methods, 51–53 have been developed and successfully applied to numerous biological systems. This encourages us to understand the robustness of Treg development through theoretical and computational exploration of the underlying regulatory network. Here, we construct a minimal transcriptional network, consisting of FoxP3, an up-regulated cofactor and an down-regulated cofactor, to mimic the complex controlling program of Treg development. We analyze the features of the development program by exploring the model’s phase diagram and illustrating its dynamic behaviors in different phases. By using the Gillespie simulation 54 as an in silico experimental method and employing the barrier heights of a potential landscape as a metric, we further show that the maintenance of multistability of the controlling network requires the orchestration of its positive and negative feedback loops and that the network


Page 4 of 32

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


robustness is combinatorially controlled by the strengths of the both loops. In contrast, the increase of system volume monotonically contributes to the network’s robustness. We further show that the model dynamics shifts in agreement with experimental findings when all of the FoxP3DNA dissociation constants are tuned proportionally. We conclude by summarizing our findings and discussing future developments of this computational framework.

Methods Model derivation and bifurcation analysis. The mathematical model (Eq. 1) was derived from the biochemical reactions involved in our minimal Treg network (Fig. 1). As shown in Tab. 1, there are mainly two classes of reactions: (1) formations and dissociations of FoxP3-cofactor and FoxP3-cofactor-promoter complexes; and (2) productions and degradations of reactant molecules, including FoxP3 and cofactors. In order to obtain the ordinary differential equation model, we have adopted a quasi-equilibrium assumption for the class (1) reactions, i.e., complex formation and dissociation occur at a much shorter time scale compared with protein production and degradation. To illustrate the dynamics of this model, MatCont, a graphical Matlab dynamics analysis package, 55 was used to perform bifurcation study for the mathematical model.

Stochastic simulations. The model (Eq. 1) offers a simple description of the dynamic and kinetic behaviors of the Treg development program but fails to offer deeper insights due to its deterministic nature. To elucidate the molecular fluctuations of the network and their contribution to system behavior, stochastic descriptions, such as Gillespie simulation and later developments, 56 can be employed. In principle, such a computational analysis can be implemented by directly simulating the biochemical reactions listed in Tab. 1. However, due to its low simulation efficiency, adoption of the original Gillespie simulation can be challenging for systems with multiple timescales, large molecular numbers, and many reaction channels. To resolve this issue without losing the stochastic nature of the system, we used an approximated simulation method that has been used in many



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

similar systems. 57 The basic idea of this approximated method is to take the overall production and degradation of a molecular species in a biochemical system as two individual reactions with their reaction rates derived from the corresponding deterministic model, instead of simulating every detailed biochemical reactions. In our case, we can derive such an approximated description from Eq. 1 as shown below: αu0 +αu Hu (Ft ,U,D)

0/ −−−−−−−−−−−→ U γuU

U −−→ 0/ αd0 +αd Hd (Ft ,U,D)

0/ −−−−−−−−−−−→ D γd D

D −−→ 0/ α f 0 +α f H f (Ft ,U,D)

0/ −−−−−−−−−−−→ Ft γ f Ft

Ft −−→ 0/

. This greatly reduces system complexity and the computational burden, enabling systematic exploration of the Treg controlling program.

Data acquisition and processing. The steady state distributions of molecular species (total FoxP3, up- and down-regulated cofactors) were obtained by performing a long stochastic simulation (> 109 simulation steps), recording the numbers of the species at equally spaced time points, and acquiring corresponding histograms by sorting the recorded data according to their values. The potential landscape of a molecular species was derived by taking the negative logarithm of the corresponding steady state distribution of the molecule. To reduce the roughness of the potential landscapes, particularly in the case when potential barrier height is low, the built-in Savitzky-Golay filter in Matlab was used for smoothing. Similarly, gaussian filtering was used for smoothing the two-dimensional potential landscape (Fig. 3F). The height difference between the bottom of a well and the top of the barrier of a landscape is defined as barrier height (Fig. 3I), which was adopted as a metric for quantifying the stability of the system in that state. 6 ACS Paragon Plus Environment

Page 6 of 32

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In addition to barrier height, mean residence time, the average amount of time that the system spends in a particular state, was also employed as an independent metric for the measurement of system robustness. To obtain the mean residence time of the system, the saddle point of the corresponding potential landscape was used as the exit boundary of a state, and a total of 103 repeats of escape simulations were performed to compute a statistical average.

Parameters. We chose parameter values through appropriate scaling of their correspondingly biologically relevant values available in BioNumbers database. 58 Specifically, the degradation rates of all of the molecules were set as unit (γu = γd = γ f = 1.0) and their basal and inducible production rate constants were set as 3.0 (αu0 , αd0 , α f 0 , αu , αd , α f ). For simplicity, we assumed the activation strengths of both FoxP3-cofactor-promoter complexes (FoxP3-U-promoter and PoxP3D-promoter) are the same for each of the positive feedback loops, i.e., C2u = C3u and C2 f = C3 f , and the basal expressions of naked and FoxP3-bounded promoter are equal in the negative feedback, i.e. (C0d = C1d ). In addition, the dissociation constants (M1 = 12, M2 = 1.2, M3 = 800, W1 = 59.2, W2 = 25.0), system volume (Ω = 2.0), and the Hill coefficients (ν = θ1 = θ2 = 2.0) were also used. We also take C0u , C1u , C0 f and C1 f all as 0.0, corresponding to zero basal production for the positive feedback induced production, and, C2d = C3d = 0.0, corresponding to a null production upon full inhibition in the negative feedback regulation. The above parameters constitute the default set of the parameter values throughout the paper unless specified elsewhere. In Fig. 2, the two-parameter bifurcation diagram (Panel A) was acquired by varying Cu and Cd , where Cu = C2u = C3u and Cd = C0d = C1d . In addition, the five sets of the parameters (Cu ,Cd ), L(1, 10), H(4.5, 0.5), B2 (3.5, 6.8), B3 (4.5, 9.6) and B1 (2.0, 3.0), were chosen for stochastic simulations. Panel B-E correspond to the first four parameter sets above. In Fig. 3, the potential landscapes in Panel A-E correspond to the parameter sets of L, H, B1 , B2 , and B3 in Fig. 2A. The B3 parameter set was used to generate Panel F and K. In Fig. 4, a total of 17 by 33 grids were scanned to produce the heatmap of the robustness metric, where a parameter set identical to that of Fig. 2A was used. In Fig. 5, all of the parameters are identical to those for Fig. 2E except the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

varying system volume. In Fig. 6, a one-parameter bifurcation diagram was obtained by scaling the dissociation constants (M1 , M2 , M3 ) proportionally (Panel A). Here, system volume was set as Ω = 1 and the rest parameters are identical to those for Fig. 2E. The four potential landscapes (a-d) correspond to a scaling of 10−1.0 , 10−0.3 , 100.3 and 102.0 for all of the dissociation constants respectively.

Results and discussion A minimal Treg model. Recent experiments have revealed that the signature of Treg can be defined by the differential expression profiles of 603 genes, compared with those of Tconv, among which 407 genes are up-regulated and the remaining 198 are down-regulated. 21 Moreover, these signature genes interact with each other to form multiple redundant regulatory feedback loops, including both positive and negative, for the auto-assembly and lock-in of the Treg state. We have mimicked such a regulatory interaction network by proposing a computational model that consists of two master species F and F*, representing the native FoxP3 and its modified functional form (i.e., functionally active form of FoxP3, such as FoxP3-cofactor complex), and a set of downstream cofactors Ui ’s and Di ’s, corresponding to the molecules of the up- and down-regulated signature genes whose expressions are significantly increased or decreased when a cell is at the Treg state 21 (Fig. 1A). This model successfully reproduced many experimental observations, including the robustness of the Treg state’s signature when faced with removal of part of the transcription factors and the partial acquisition of the Treg signature in the Treg-like cells. However, the complexity of the model prohibits analytical analysis of the underlying network, failing to provide mechanistic insights of the controlling network of Treg development. Moreover, the deterministic description of the model overlooks the molecular fluctuations intrinsic to the system, rendering the inability in addressing the stochasticity and functional robustness of the Treg network. To enable a systematic exploration of the controlling network of Treg development, we propose a minimal transcriptional network that consists of two master species F and F*, an up-regulated


Page 8 of 32

Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cofactor U, and a down-regulated cofactor D as shown in Fig. 1B. Here, FoxP3 transits from its native (F) to functional (F*) state, mediated by its cofactors (U and D), while its functional form autonomously promotes the production of the native state. In addition, although both forms of FoxP3 regulate the cofactor genes, the modified one has a much stronger regulation than the native form. This minimal model enables mathematical analysis but still retains the key features of the Treg controlling program, enabling us to investigate the program’s dynamic behavior and system characteristics. The corresponding biochemical reactions of the minimal network are listed in Tab. 1, from which we can derive a mathematical model (see the Methods section for details) as dU = αu0 + αu Hu (Ft ,U, D) − γuU dt dD = αd0 + αd Hd (Ft ,U, D) − γd D dt dFt = α f 0 + α f H f (Ft ,U, D) − γ f Ft dt

(1)

where Ft is total concentration of native FoxP3 (F) and its modified complex (F*), U and D are the concentrations of the up- and down-regulated cofactor genes that simultaneously serve as transcriptional activators of FoxP3 gene. αs0 , αs and γs (s = u, d, f ) are the rate constants for the basal production, inducible production and degradation of a gene respectively. The function Hs (Ft ,U, D) (s = u, d, f ) refers to the hybrid transcription kinetics that are co-regulated by Ft , U and D !υ C0,s + H{s=u,d, f } (Ft ,U, D) =

Ft θ

θ

1

2

1+( UW 1 )+( DW 2 )

!υ 1+

n

Ft θ

θ

1

2

1+( UW 1 )+( DW 2 )

C1,s M1

θ C + M2,s2 ( UW11 )υ

θ C + M3,s3 ( DW22 )υ

o (2)

n

1 M1

U θ1

Dθ2

+ M12 ( W1 )υ + M13 ( W2 )υ

o

where C0,s , C1,s , C2,s , C3,s are the folds of change for the expression level of species s (s = u, d, f ) when promoters are bound by nothing, original FoxP3, FoxP3-U complex, and FoxP3-D complex respectively. The dissociation constants W1,2 and M1,2,3 and the Hill coefficients θ1,2 and υ are detailed in Tab. 1.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Phase diagram and stochastic network dynamics. The simplicity of the minimal model enables us to systematically explore its dynamic features by performing bifurcation analysis. As depicted in the phase diagram (Fig. 2A), the Treg network may exhibit distinct dynamic behaviors depending on its parameters: The system (total FoxP3 level here) could be monostably low despite its initial conditions (upper left region of the diagram), monostably high (lower right region), or bistable (shaded upper right region). Bistability means that total FoxP3 expression level can be either stably high or stably low, corresponding to the two states of T cells - Treg and Tconv. To illustrate the dynamic manifold of the Treg network, we constructed a stochastic model by adopting the Gillespie algorithm as our simulation core (Detailed in Methods Section). Different from the deterministic description (Eq. 1) where noise is eliminated, this stochastic model intrinsically takes molecular fluctuations into account and thus enables us to uncover more realistic system dynamics. Fig. 2B-E show representative time evolution trajectories of the total FoxP3 level. It is clear that the Treg networks may remain stably low (Fig. 2B) or stably high (Fig. 2C), or exhibit bistable behaviors (Fig. 2D and 2E), regardless of the intrinsic fluctuations of molecular numbers. In addition, Fig. 2D and 2E also show that T cells can switch spontaneously from one state to the other without external perturbations. The transition of cellular states is purely driven by the fluctuations of molecules. Comparison of Fig. 2D and 2E further suggests that the occurrence of state transition events within a given amount of time, i.e., the rate of state transition, could be dramatically different, even though the parameters of both cases are within the bistability regime. As the low and high states of the model correspond to the Tconv and Treg states of T cells, this rate of transition reflects how stably the system maintains its state, which encourages us to quantify the functional robustness of the Treg controlling network by characterizing its state transition frequency.

Potential landscape and robustness metric. To enable systematic investigation of the stability of the Treg controlling network, we introduce the energy landscape concept for our model. An energy landscape is a mapping of all possible conformations of a molecular entity, or the spatial


Page 10 of 32

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


positions of interacting molecules in a system, and their corresponding energy levels, onto a twoor three-dimensional Cartesian coordinate system. 59 Such a mapping has been widely used in several biological problems, such as protein folding. 60,61 As most cellular networks, such as those for signal transduction, molecular interactions, and cellular development including Treg development studying here, are out of equilibrium, a system’s landscape can not be directly obtained from its hamiltonian or energy function. Instead, it needs to be derived from the steady state probability distribution of the system. 62–72 Such an effective potential landscape approach intrinsically encodes the dynamic property (e.g., system stability) of a system into its landscape profile. Mathematically, this can be implemented by taking the negative logarithm of the corresponding steady state probability distribution, i.e., Φ(x) = − log P(x), where P(x) is the steady state probability distribution of molecules x (any of total FoxP3, U, and D or their combinations). The corresponding steady state distributions of molecular species can be acquired by sorting data recorded from a long run of stochastic simulation (See the Methods section). Fig. 3A-E show the potential landscapes of the Treg network under different parameter regimes, corresponding to the dots L, H, B1 , B2 , and B3 in the phase diagram (Fig. 2A). Although the actual potential landscape is high-dimensional (total FoxP3, up- and down-regulated cofactors), we have projected it onto the total FoxP3 axis for a simpler presentation. We chose total FoxP3 as the variable of the landscape because of its central role to the functionality of Treg development. To examine if such a degenerated presentation indeed represents a correct multistability picture of the system, we also plotted the network’s potential landscape in the space of the two cofactors (U and D) for the case of Fig. 3E. As shown in Fig. 3F, the two-dimensional landscape indeed has two wells corresponding to two stable states, consistent with the double-well characteristics of Fig. 3E. The corresponding landscapes can be derived following the procedures as schematically illustrated in Fig. 3G-H. As Kramer’s theory states that the transition rate of a system from one state to another is exponentially proportional to the negative value of the barrier height of the system’s potential landscape 73 (Fig. 3I), the landscape view of the Treg network offers a simple solution for quantitative



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

measurement of the system’s robustness. Specifically, we propose to use the geometric mean of √ the two barrier heights, i.e., H1 H2 , as a robustness metric, where H1 and H2 are the heights of the two wells. We chose the geometric mean, not any individual, of the two heights for quantification is due to the fact that a functionally ‘correct’ Treg network needs to have two states, corresponding to Tconv and Treg, and the stabilities of the network in both states are important to the system’s robustness. To illustrate this idea, we analyzed the barrier heights for the above five sample cases in Fig. 3J, where the robustness metric increases with averaged barrier heights, matching well with reduced transition frequency illustrated in the simulation results in Fig. 2. To √ further validate if H1 H2 is a good robustness metric, we evaluated the mean residence time of √ the system in each state (τ1 and τ2 ) (Fig. 3K) and used their geometric mean ( τ1 τ2 ) to quantify system stability. Fig. 3L indeed shows a monotonic increase of mean residence time with barrier height, qualitatively consistent with Fig. 3J. It is worthy to note that there are other approaches for robustness quantification, for example, the least dissipation cost required for budding yeast to maintain its cell cycle, could be an alternative robustness metric. 74 In addition, rigorously speaking, barrier heights used for robustness measurement need to be derived from the actual high-dimensional potential landscape. However, the one-dimensional landscape based metric proposed here indeed allows for a simple but yet quantitative characterization of network robustness.

Feedback strengths to network robustness. Recent experimental studies have shown that feedback loops, both positive and negative, lie in the center of the architecture of the Treg controlling network . 75–77 To elucidate their roles in contributing to the network’s functional robustness, we conducted a computational survey by performing stochastic simulations for different combinations of feedback strengths and evaluating corresponding system robustness by using the robustness metric as a measurement. Fig. 4 shows the network’s robustness as a function of the loops’ strength. It is clear that the system can lose its bistability with monotonic increase or decrease of the strength of one feedback loop while keeping the other constant (the broken lines 1 and 2). In contrast,


Page 12 of 32

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


coordinated tuning of the feedback strengths (the broken line 3) may lead to increased network robustness. This finding can be understood by considering limiting cases: With a given strength of the negative feedback loop, the system may become monostable (monostably high) when the positive feedback loop is infinitely strong; conversely, a system with a given positive feedback strength will be shifted to be monostably low by an infinite strength of negative feedback. As well documented in literature, 78–80 positive feedback is essential and self-sufficient for the generation of a system’s multistability, and is also important for producing irreversibility of a system’s state. In contrast, single negative feedback does not result in multistability (double negative feedback does), but it contributes to system stability by limiting the range over which the concentrations of network components fluctuate. 32,80 This result suggests that the functional robustness of the Treg developmental program is not determined solely by a single feedback loop, but rather the orchestration of the both, which indicates an architectural requirement for the accomplishment of stable state maintenance during Treg development. From a network designing perspective, such a requirement sets an architectural constraint over the network for desired system functionality.

Role of system volume to network robustness. Through the combination of stochastic simulation and robustness measurement, we also investigated the role of system volume in the Treg network’s robustness. As illustrated in Fig. 5A, the robustness of the Treg network increases monotonically with system volume. In concert with that, the barrier heights of the corresponding potential landscapes increase monotonically as well (Ω = 0.5, 1.0, 1.5, 2.0, 2.5) (Fig. 5B). This can be interpreted with the following argument: The number variance of a molecular species is proportional to the root square of its mean, which suggests that the relative amplitude of fluctuations (variance divided by its mean) is inversely proportional to the root square of the mean molecular number. On the other hand, a larger system corresponds to a larger absolute number of molecules because the molecular number of a species equals to its concentration multiplied by system volume. As a result, fluctuations of molecules are inversely proportional to the root square of system volume. This relation suggests that an increased system volume results in a



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

smaller molecular fluctuation, which in turn leads to a rarer transition and hence increased system robustness.

FoxP3-DNA binding affinity to network robustness. Supporting the role of FoxP3 as the core regulator of the Treg network for the state development and maintenance, recent loss-of-function studies have shown that FoxP3 mutations can lead to the defective development of functional Treg cells and that such a defection can further result in a fatal autoimmune disorder in both mice and humans. 9,10 To find possible constraints of FoxP3 for appropriate network functioning, we employed the same set of simulation programs and measurement protocols to study the relationship between the functional robustness of the Treg network and the kinetic features of FoxP3, particularly the binding affinity of FoxP3 with the promoters that it regulates. Specifically, we mimicked experimental FoxP3 mutations by tuning all of the dissociation constants (M1 , M2 , M3 ) proportionally with a scaling factor defined as M. Fig. 6A depicts the steady state value of total FoxP3 with respect to the scaling factor, obtained through the bifurcation analysis of the deterministic model Eq. 1. Clearly, the system varies from a single stable state phase (high total FoxP3, corresponding to the Treg state) to a bistable phase, corresponding to the coexistence of Tconv and Treg states, as the scaling factor increases. The system can further turn into a monostably low state (low total FoxP3, corresponding to the Tconv state) when the scaling factor is large. The alteration of the system dynamics can also be illustrated with the potential landscape language: as shown in Fig. 6A, the four plots at the bottom indicate distinct landscapes for different scaling factors. In addition to mimicking the loss-of-function experiment, computational probing of the FoxP3DNA binding affinity may provide insights for signal-induced transition of CD4+ T cells from the Tconv to the Treg states. Several experiments have shown that both the FoxP3 locus and the loci to which FoxP3 binds are subject to epigenetic modifications (i.e., methylation and demethylation), 81–83 indicating that appropriate external signals might be able to alter FoxP3-DNA binding affinity and hence the system’s potential landscape profile. If the commitment to Treg state is


Page 14 of 32

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


associated with such signals, the robustness of cellular states could be greatly reinforced by the accompanying alteration of potential landscapes if the system has a landscape profile similar to Fig. 6A(c) in the Tconv state but is shifted to have a profile like Fig. 6A(b) during the induction. Such a hypothesis is in line with Hori’s speculation that stable FoxP3 expression is ensured by (1) demethylation of a Treg-specific region independent of FoxP3 and (2) a positive feedback mechanism through which FoxP3 subsequently stabilizes its own transcription. 84 Although having not been directly demonstrated by experimental reports, this hypothesis indeed raises an interesting possible mechanism for robust control over Treg development and needs experimental validation in the near future.

Conclusion In this paper, we have constructed a minimal transcriptional network that consists of FoxP3, one up- and one down-regulated cofactors, to mimic the Treg controlling program and used this model to systematically examine the functional robustness of Treg development by combining dynamics analysis with stochastic simulation. We found that feedback strengths, system volume, and FoxP3DNA binding affinity all contribute to the network’s robustness despite their difference in specific impacts. This work provides a system-level view of and a quantitative understanding of the Treg developmental program, and also offers new insights into the dynamics and robustness features of regulatory T cells. In addition, recent efforts have shown that the transfer of regulatory T cells can be a promising therapy that holds tremendous application potentials, including specific suppression of autoimmunity, treatment of graft versus-host disease, and improvement of acceptance of transplanted organs. 85,86 Therefore, our study of the minimal network is also valuable for the development of future therapies. The quantitative and systems understanding of the Treg controlling network has just started with the advance of experimental efforts, and many important aspects need to be investigated in



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

near future. For instance, the actual Treg network possesses multiple redundant feedback loops that were not included in our minimal model for simplicity; however, such an architectural redundancy can be crucial for the functional robustness of the developmental program. In addition, our current study has only explored the roles of interaction strengths to Treg development, but the architectural topology of the controlling network is equally important to system behavior. Moreover, the roles of triggering factors that induce cellular state transitions need to be explored as well. Nevertheless, the minimal model presented here and associated analysis method may serve as a versatile computational platform for quantitative understanding of Treg development, offering both fundamental and therapeutic implications.

Acknowledgement This work was supported in part by American Heart Association under Grant No. 12SDG12090025. We thank Andrew Blanchard for commenting and editing the manuscript.

References (1) Belkaid, Y.; Rouse, B. T. Natural Regulatory T Cells in Infectious Disease. Nat. Immunol. 2005, 6, 353–360. (2) Sakaguchi, S.; Yamaguchi, T.; Nomura, T.; Ono, M. Regulatory T Cells and Immune Tolerance. Cell 2008, 133, 775–787. (3) Miyara, M.; Sakaguchi, S. Natural Regulatory T Cells: Mechanisms of Suppression. Trends Mol. Med. 2007, 13, 108–116. (4) Josefowicz, S. Z.; Lu, L.-F.; Rudensky, A. Y. Regulatory T Cells: Mechanisms of Differentiation and Function. Annu. Rev. Immunol. 2012, 30, 531–564. (5) Fontenot, J. D.; Rudensky, A. Y. A Well Adapted Regulatory Contrivance: Regulatory T Cell


Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Development and the Forkhead Family Transcription Factor Foxp3. Nat. Immunol. 2005, 6, 331–337. (6) Ramsdell, F. Foxp3 and Natural Regulatory T Cells: Key to a Cell Lineage? Immunity 2003, 19, 165–168. (7) Zheng, Y.; Rudensky, A. Y. Foxp3 in Control of the Regulatory T Cell Lineage. Nat. Immunol. 2007, 8, 457–462. (8) Hori, S.; Nomura, T.; Sakaguchi, S. Control of Regulatory T Cell Development by the Transcription Factor Foxp3. Sci. Signal. 2003, 299, 1057–1061. (9) Fontenot, J. D.; Gavin, M. A.; Rudensky, A. Y. Foxp3 Programs the Development and Function of CD4+ CD25+ Regulatory T Cells. Nat. Immunol. 2003, 4, 330–336. (10) Khattri, R.; Cox, T.; Yasayko, S.-A.; Ramsdell, F. An Essential Role for Scurfin in CD4+ CD25+ T Regulatory Cells. Nat. Immunol. 2003, 4, 337–342. (11) Williams, L. M.; Rudensky, A. Y. Maintenance of the Foxp3-Dependent Developmental Program in Mature Regulatory T Cells Requires Continued Expression of Foxp3. Nat. Immunol. 2007, 8, 277–284. (12) Zheng, Y.; Josefowicz, S. Z.; Kas, A.; Chu, T.-T.; Gavin, M. A.; Rudensky, A. Y. GenomeWide Analysis of Foxp3 Target Genes in Developing and Mature Regulatory T Cells. Nature (London) 2007, 445, 936–940. (13) Marson, A.; Kretschmer, K.; Frampton, G. M.; Jacobsen, E. S.; Polansky, J. K.; MacIsaac, K. D.; Levine, S. S.; Fraenkel, E.; Von Boehmer, H.; Young, R. A. Foxp3 Occupancy and Regulation of Key Target Genes During T-Cell Stimulation. Nature (London) 2007, 445, 931–935. (14) Gavin, M. A.; Rasmussen, J. P.; Fontenot, J. D.; Vasta, V.; Manganiello, V. C.; Beavo, J. A.;



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Rudensky, A. Y. Foxp3-Dependent Programme of Regulatory T-Cell Differentiation. Nature (London) 2007, 445, 771–775. (15) Lin, W.; Haribhai, D.; Relland, L. M.; Truong, N.; Carlson, M. R.; Williams, C. B.; Chatila, T. A. Regulatory T Cell Development in the Absence of Functional Foxp3. Nat. Immunol. 2007, 8, 359–368. (16) Pfoertner, S.; Jeron, A.; Probst-Kepper, M.; Guzman, C. A.; Hansen, W.; Westendorf, A. M.; Toepfer, T.; Schrader, A. J.; Franzke, A.; Buer, J.; et al. Signatures of Human Regulatory T Cells: An Encounter with Old Friends and New Players. Genome Biol. 2006, 7, R54. (17) Birzele, F.; Fauti, T.; Stahl, H.; Lenter, M. C.; Simon, E.; Knebel, D.; Weith, A.; Hildebrandt, T.; Mennerich, D. Next-Generation Insights into Regulatory T Cells: Expression Profiling and FoxP3 Occupancy in Human. Nucleic Acids Res. 2011, 39, 7946–7960. (18) Hill, J. A.; Feuerer, M.; Tash, K.; Haxhinasto, S.; Perez, J.; Melamed, R.; Mathis, D.; Benoist, C. Foxp3 Transcription-Factor-Dependent and -Independent Regulation of the Regulatory T Cell Transcriptional Signature. Immunity 2007, 27, 786–800. (19) Delgoffe, G. M.; Bettini, M. L.; Vignali, D. A. Identity Crisis: It’s Not Just Foxp3 Anymore. Immunity 2012, 37, 759–761. (20) Rudra, D.; Chaudhry, A.; Niec, R. E.; Arvey, A.; Samstein, R. M.; Leslie, C.; Shaffer, S. A.; Goodlett, D. R.; Rudensky, A. Y. Transcription Factor Foxp3 and Its Protein Partners Form a Complex Regulatory Network. Nat. Immunol. 2012, 13, 1010–1019. (21) Fu, W.; Ergun, A.; Lu, T.; Hill, J. A.; Haxhinasto, S.; Fassett, M. S.; Gazit, R.; Adoro, S.; Glimcher, L.; Chan, S.; et al. A Multiply Redundant Genetic Switch ‘Locks in’ the Transcriptional Signature of Regulatory T Cells. Nat. Immunol. 2012, 13, 972–980. (22) Hori, S. The Foxp3 Interactome: A Network Perspective of Treg Cells. Nat. Immunol. 2012, 13, 943–945. 18 ACS Paragon Plus Environment

Page 18 of 32

Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(23) Rubtsov, Y. P.; Niec, R. E.; Josefowicz, S.; Li, L.; Darce, J.; Mathis, D.; Benoist, C.; Rudensky, A. Y. Stability of the Regulatory T Cell Lineage in Vivo. Science 2010, 329, 1667–1671. (24) Bailey-Bucktrout, S. L.; Bluestone, J. A. Regulatory T Cells: Stability Revisited. Trends Immunol. 2011, 32, 301–306. (25) Hamann, A. Regulatory T Cells Stay on Course. Immunity 2012, 36, 161–163. (26) Miyao, T.; Floess, S.; Setoguchi, R.; Luche, H.; Fehling, H. J.; Waldmann, H.; Huehn, J.; Hori, S. Plasticity of Foxp3 T Cells Reflects Promiscuous Foxp3 Expression in Conventional T Cells but Not Reprogramming of Regulatory T Cells. Immunity 2012, 36, 262–275. (27) Mendoza, L.; Pardo, F. A Robust Model to Describe the Differentiation of T-Helper Cells. Theor. Biosci. 2010, 129, 283–293. (28) Hong, T.; Xing, J.; Li, L.; Tyson, J. J. A Mathematical Model for the Reciprocal Differentiation of T Helper 17 Cells and Induced Regulatory T Cells. PLoS Comput. Biol. 2011, 7, e1002122. (29) Hong, T.; Xing, J.; Li, L.; Tyson, J. A Simple Theoretical Framework for Understanding Heterogeneous Differentiation of CD4+ T Cells. BMC Syst. Biol. 2012, 6, 66. (30) Roeder, I.; Glauche, I. Towards an Understanding of Lineage Specification in Hematopoietic Stem Cells: A Mathematical Model for the Interaction of Transcription Factors GATA-1 and PU.1. J. Theor. Biol. 2006, 241, 852–865. (31) Kitano, H. Biological Robustness. Nat. Rev. Genet. 2004, 5, 826–837. (32) Stelling, J.; Sauer, U.; Szallasi, Z.; Doyle, F. J.; Doyle, J. Robustness of Cellular Functions. Cell 2004, 118, 675–685. (33) Barkai, N.; Shilo, B.-Z. Variability and Robustness in Biomolecular Systems. Mol. Cell. 2007, 28, 755–760. 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(34) Li, F.; Long, T.; Lu, Y.; Ouyang, Q.; Tang, C. The Yeast Cell-Cycle Network is Robustly Designed. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 4781–4786. (35) Little, J. W.; Shepley, D. P.; Wert, D. W. Robustness of a Gene Regulatory Circuit. EMBO J. 1999, 18, 4299–4307. (36) Shah, N. A.; Sarkar, C. A. Robust Network Topologies for Generating Switch-Like Cellular Responses. PLoS Comput. Biol. 2011, 7, e1002085. (37) Venturelli, O. S.; El-Samad, H.; Murray, R. M. Synergistic Dual Positive Feedback Loops Established by Molecular Sequestration Generate Robust Bimodal Response. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, E3324–E3333. (38) Kærn, M.; Elston, T. C.; Blake, W. J.; Collins, J. J. Stochasticity in Gene Expression: From Theories to Phenotypes. Nat. Rev. Genet. 2005, 6, 451–464. (39) Raj, A.; van Oudenaarden, A. Nature, Nurture, or Chance: Stochastic Gene Expression and Its Consequences. Cell 2008, 135, 216–226. (40) Balázsi, G.; van Oudenaarden, A.; Collins, J. J. Cellular Decision Making and Biological Noise: From Microbes to Mammals. Cell 2011, 144, 910–925. (41) Swain, P. S.; Elowitz, M. B.; Siggia, E. D. Intrinsic and Extrinsic Contributions to Stochasticity in Gene Expression. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12795–12800. (42) Lu, T.; Ferry, M.; Weiss, R.; Hasty, J. A Molecular Noise Generator. Phys. Biol. 2008, 5, 036006. (43) Rao, C. V.; Wolf, D. M.; Arkin, A. P. Control, Exploitation and Tolerance of Intracellular Noise. Nature (London) 2002, 420, 231–237. (44) Zhuravlev, P. I.; Papoian, G. A. Molecular Noise of Capping Protein Binding Induces Macroscopic Instability in Filopodial Dynamics. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 11570– 11575. 20 ACS Paragon Plus Environment

Page 20 of 32

Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(45) Raser, J. M.; O’Shea, E. K. Noise in Gene Expression: Origins, Consequences, and Control. Science 2005, 309, 2010–2013. (46) Sasai, M.; Wolynes, P. G. Stochastic Gene Expression as a Many-Body Problem. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 2374–2379. (47) Aurell, E.; Sneppen, K. Epigenetics as a First Exit Problem. Phys. Rev. Lett. 2002, 88, 48101. (48) Wang, J.; Zhang, K.; Wang, E. Kinetic Paths, Time Scale, and Underlying Landscapes: A Path Integral Framework to Study Global Natures of Nonequilibrium Systems and Networks. J. Chem. Phys. 2010, 133, 125103. (49) Assaf, M.; Roberts, E.; Luthey-Schulten, Z. Determining the Stability of Genetic Switches: Explicitly Accounting for mRNA Noise. Phys. Rev. Lett. 2011, 106, 248102. (50) Zong, C.; Lu, T.; Shen, T.; Wolynes, P. G. Nonequilibrium Self-Assembly of Linear Fibers: Microscopic Treatment of Growth, Decay, Catastrophe and Rescue. Phys. Biol. 2006, 3, 83– 92. (51) Tian, T.; Burrage, K. Stochastic Models for Regulatory Networks of the Genetic Toggle Switch. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 8372–8377. (52) Shen, T.; Hamelberg, D. A Statistical Analysis of the Precision of Reweighting-Based Simulations. J. Chem. Phys. 2008, 129, 034103. (53) Lu, T.; Volfson, D.; Tsimring, L.; Hasty, J. Cellular Growth and Division in the Gillespie Algorithm. Syst. Biol. 2004, 1, 121–128. (54) Gillespie, D. T. Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys. Chem. 1977, 81, 2340–2361. (55) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A. MATCONT: A MATLAB Package for Numerical Bifurcation Analysis of ODEs. ACM TOMS 2003, 29, 141–164. 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(56) Gillespie, D. T. Stochastic Simulation of Chemical Kinetics. Annu. Rev. Phys. Chem. 2007, 58, 35–55. (57) Süel, G. M.; Garcia-Ojalvo, J.; Liberman, L. M.; Elowitz, M. B. An Excitable Gene Regulatory Circuit Induces Transient Cellular Differentiation. Nature (London) 2006, 440, 545–550. (58) Milo, R.; Jorgensen, P.; Moran, U.; Weber, G.; Springer, M. BioNumbers–The Database of Key Numbers in Molecular and Cell Biology. Nucleic Acids Res. 2010, 38, D750–D753. (59) Wales, D. Energy Landscapes: Applications to Clusters, Biomolecules and Glasses; Cambridge University Press, 2004. (60) Zheng, W.; Schafer, N. P.; Davtyan, A.; Papoian, G. A.; Wolynes, P. G. Predictive Energy Landscapes for Protein–Protein Association. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 19244– 19249. (61) Onuchic, J. N.; Luthey-Schulten, Z.; Wolynes, P. G. THEORY OF PROTEIN FOLDING: The Energy Landscape Perspective. Annu. Rev. Phys. Chem. 1997, 48, 545–600. (62) Lapidus, S.; Han, B.; Wang, J. Intrinsic Noise, Dissipation Cost, and Robustness of Cellular Networks: The Underlying Energy Landscape of MAPK Signal Transduction. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 6039–6044. (63) Wang, J.; Xu, L.; Wang, E.; Huang, S. The Potential Landscape of Genetic Circuits Imposes the Arrow of Time in Stem Cell Differentiation. Biophys. J. 2010, 99, 29–39. (64) Wang, J.; Li, C.; Wang, E. Potential and Flux Landscapes Quantify the Stability and Robustness of Budding Yeast Cell Cycle Network. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 8195–8200. (65) Wang, J.; Zhang, K.; Xu, L.; Wang, E. Quantifying the Waddington Landscape and Biological Paths for Development and Differentiation. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 8257– 8262. 22 ACS Paragon Plus Environment

Page 22 of 32

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(66) Zhang, F.; Xu, L.; Zhang, K.; Wang, E.; Wang, J. The Potential and Flux Landscape Theory of Evolution. J. Chem. Phys. 2012, 137, 065102. (67) Cao, Y.; Lu, H.-M.; Liang, J. Probability Landscape of Heritable and Robust Epigenetic State of Lysogeny in Phage Lambda. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 18445–18450. (68) Strasser, M.; Theis, F. J.; Marr, C. Stability and Multiattractor Dynamics of a Toggle Switch Based on a Two-Stage Model of Stochastic Gene Expression. Biophys. J. 2012, 102, 19–29. (69) Kim, K.-Y.; Wang, J. Potential Energy Landscape and Robustness of a Gene Regulatory Network: Toggle Switch. PLoS Comput. Biol. 2007, 3, e60. (70) Wang, J.; Xu, L.; Wang, E. Potential Landscape and Flux Framework of Nonequilibrium Networks: Robustness, Dissipation, and Coherence of Biochemical Oscillations. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 12271–12276. (71) Sisan, D. R.; Halter, M.; Hubbard, J. B.; Plant, A. L. Predicting Rates of Cell State Change Caused by Stochastic Fluctuations Using a Data-Driven Landscape Model. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 19262–19267. (72) Zhou, J. X.; Aliyu, M.; Aurell, E.; Huang, S. Quasi-Potential Landscape in Complex MultiStable Systems. J. R. Soc. Interface 2012, 9, 3539–3553. (73) Mel’Nikov, V. The Kramers Problem: Fifty Years of Development. Phys. Rep. 1991, 209, 1–71. (74) Han, B.; Wang, J. Least Dissipation Cost as a Design Principle for Robustness and Function of Cellular Networks. Phys. Rev. E 2008, 77, 031922. (75) Qi, H.; Blanchard, A.; Lu, T. Engineered Genetic Information Processing Circuits. WIREs Syst. Biol. Med. 2013, 5, 273–287. (76) Becskei, A.; Serrano, L. Engineering Stability in Gene Networks by Autoregulation. Nature (London) 2000, 405, 590–593. 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(77) Acar, M.; Becskei, A.; van Oudenaarden, A. Enhancement of Cellular Memory by Reducing Stochastic Transitions. Nature (London) 2005, 435, 228–232. (78) Xiong, W.; Ferrell Jr, J. E. A Positive-Feedback-Based Bistable ‘Memory Module’ that Governs a Cell Fate Decision. Nature (London) 2003, 426, 460–465. (79) Ferrell Jr, J. E. Self-Perpetuating States in Signal Transduction: Positive Feedback, DoubleNegative Feedback and Bistability. Curr. Opin. Chem. Biol. 2002, 6, 140–148. (80) Tyson, J. J.; Chen, K. C.; Novak, B. Sniffers, Buzzers, Boggles and Blinkers: Dynamics of Regulatory and Signaling Pathways in the Cell. Curr. Opin. Cell Biol. 2003, 15, 221–231. (81) Floess, S.; Freyer, J.; Siewert, C.; Baron, U.; Olek, S.; Polansky, J.; Schlawe, K.; Chang, H.D.; Bopp, T.; Schmitt, E.; et al. Epigenetic Control of the Foxp3 Locus in Regulatory T Cells. PLoS Biol. 2007, 5, e38. (82) Chen, C.; Rowell, E. A.; Thomas, R. M.; Hancock, W. W.; Wells, A. D. Transcriptional Regulation by Foxp3 is Associated with Direct Promoter Occupancy and Modulation of Histone Acetylation. J. Biol. Chem. 2006, 281, 36828–36834. (83) Huehn, J.; Polansky, J. K.; Hamann, A. Epigenetic Control of FOXP3 Expression: The Key to a Stable Regulatory T-Cell Lineage? Nat. Rev. Immunol. 2009, 9, 83–89. (84) Hori, S. Developmental Plasticity of Foxp3+ Regulatory T Cells. Curr. Opin. Immunol. 2010, 22, 575–582. (85) Xu, W.; Lan, Q.; Chen, M.; Chen, H.; Zhu, N.; Zhou, X.; Wang, J.; Fan, H.; Yan, C.-S.; Kuang, J.-L.; et al. Adoptive Transfer of Induced-Treg Cells Effectively Attenuates Murine Airway Allergic Inflammation. PLoS ONE 2012, 7, e40314. (86) Komanduri, K. V.; Champlin, R. E. Can Treg Therapy Prevent GVHD? Blood 2011, 117, 751–752.


Page 24 of 32

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Biochemical reactions underlying the transcriptional controlling network of Treg development. Here, F, O, and OFν refer to FoxP3, promoter, and FoxP3-promoter complex; Xi ’s (i = 1, 2) are cofactors with X1 and X2 corresponding to the proteins of the up- (U) and down-regulated (D) genes respectively; and FXiθi and O(FXiθi )ν are FoxP3-cofactor complex and FoxP3-cofactorpromoter complex correspondingly. Class 1 reaction

Disso. const.

Explanation

[F] + θi [Xi ] [FXiθi ] [O] + ν[F] [OFν ] [OFν ] + νθi [Xi ] [O(FXiθi )ν ] [O] + ν[FXiθi ] [O(FXiθi )ν ]

Wi M1 Wi Mi+1

FoxP3 reversibly binds to signature species to form a complex Promoter and FoxP3 reversibly form a complex Promoter-FoxP3 complex and protein Xi reversibly form a complex Promoter and FoxP3-protein complex reversibly form a complex

Class 2 reaction

Rate const.

Explanation

[O] → [O] + [F] [OFν ] → [OFν ] + [F] [O(FXiθi )ν ] → [O(FXiθi )ν ] + [F]

α f C0 f α f C1 f α f C(i+1) f αf0 γf αx j C0x j αx j C1x j αx j C(i+1)x j αx j 0 γx j

Naked promoter produces FoxP3 Promoter bounded with FoxP3 produces FoxP3 Promoter bounded with FoxP3-protein Xi complex produces FoxP3.

0/ → [F] [F] → 0/ [O] → [O] + [X j ] [OFν ] → [OFν ] + [X j ] [O(FXiθi )ν ] → [O(FXiθi )ν ] + [X j ] 0/ → [X j ] [X j ] → 0/

basal-level production of FoxP3 FoxP3 degradation Naked promoter produces X j Promoter bounded with FoxP3 produces protein X j Promoter bounded with FoxP3-protein Xi complex produces protein X j basal-level production of protein X j Protein X j degradation



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

A

B

U1

F

Page 26 of 32

U2 U3

F

U

F*

D

U4

F*

D3 D2 D1

Figure 1: A mathematical model of regulatory T cell development. (A) A comprehensive model that consists of FoxP3 and multiple regulatory cofactors (Ui ’s and Di ’s) which collectively form multiple regulatory feedback loops. (B) A minimal transcriptional controlling network consisting of FoxP3, an up-regulated cofactor U, and a down-regulated cofactor D.


12

L

B3

8

60

B2

6

B1 H 1

2

1.0

D 90

2

0

30

0 0

3

Cu

4

5

6

FoxP3 total

4

0

90

Time

2.0

x105

90

60

E

30

1.0

Time

2.0

x105

3.0

30

0 0

3.0

60

0 0

C

FoxP3 total

10

B

FoxP3 total

A

Cd

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FoxP3 total

Page 27 of 32

1.0

90

Time

2.0

x105

3.0

60

30

0 0

1.0

Time

2.0

x105

3.0

Figure 2: Phase diagram and stochastic network dynamics. (A) Phase diagram. The Treg controlling network exhibits distinct behaviors in different parameter regimes, including monostably low FoxP3 (Tconv state), monostably high FoxP3 (Treg state), and bistable FoxP3 levels, corresponding to the upper left, lower right, and shaded area of the diagram. (B-E) Sample trajectories of the stochastic FoxP3 dynamics in different parameter regimes. Panel B to E correspond to the parameter sets at the dots L, H, B2 and B3 in Panel A respectively.



5 0 0

10

20

30

Foxp3 total (Ft )

40

10

5 0 0

10

20

30

Foxp3 total (Ft )

40

10

5 0 0

E

15

10

20

30

Foxp3 total (Ft )

40

15

Potential landscape

10

D

15

Potential landscape

Potential landscape

10

C

15

Potential landscape

B

15

Potential landscape

A

10

5 0 0

10

20

30

Foxp3 total (Ft )

40

5 0 0

10

20

30

Foxp3 total (Ft )

40

6

20 40 Down-regulated TF [D]

60

Probability

Time

J+

2

H1 H2 H1*H2

9

Probability

Barrier height

8

K

6 3 0

B

C

D

E

1

H1

H2 Foxp3 total

Potential landscape

L

−4

1

0

Potential landscape

Potential landscape

Histogram exponential fit

0 A

x10

I

2 3 4 x10 4 Residence time

6 Mean Residence Time

0 0

10

Probability

Time

H FoxP3 total

14

15

FoxP3 total

18

30

FoxP3 total

22

FoxP3 total

45

FoxP3 total

FoxP3 total

G F

Up-regulated TF [U]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 32

x10 3

τ τ τ1*τ 1 2

4

2

2

0

C

D

E

Figure 3: Potential landscapes of the minimal network and its application to robustness measurement. (A-E) Potential landscape of the network as a function of total FoxP3. Panel A to E correspond to the parameter sets at the dots L, H, B1 , B2 and B3 in Fig. 2A respectively. (F) The potential landscape of the network as a function of up- and down-regulated transcriptional factors, with its parameter set identical to that of Panel E. (G-H) Schematic illustration of the procedure for deriving the potential landscape. System robustness (e.g., Panel G vs. H) is intrinsically encoded in its corresponding landscape profile. (I) A schematic plot showing the barrier heights (H1 and H2 ) of a double-well landscape. (J) Robustness measurement by using the geometric mean of barrier heights as a metric. (K) Statistics of the residence time of Treg system in the low state (Tconv state) for the case corresponding to Panel E. (L) Robustness measurement by using the geometric mean of averaged residence time as a metric.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


8

1 3

6

Cd

Page 29 of 32

6

4

4

2

2

2 0 0

1

2

Cu

3

4

0

Figure 4: Network robustness with respect to feedback strengths. The orchestration of the strengths of both the negative and positive feedback loops (Cu and Cd ) determines the network robustness. Here, color indicates the degree of robustness.



12

Barrier Height

A

8

4

H1 H2 H1*H2

0 0.5

B

15

Potential landscape

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1.0 1.5 2.0 System Volume (Ω)

2.5

Ω=2.5 Ω=2.0 Ω=1.5 Ω=1.0 Ω=0.5

10

5

0

0

10

20 30 Foxp3 total (Ft )

40

Figure 5: Network robustness with respect to system volume. (A) Monotonic increase of the network’s robustness as a function of system volume. (B) Representative potential landscapes with different system volumes.


Page 30 of 32

A

30

a

b

c

B+

d

18

20

12

10

6 0 −2 10

10

−1

0

1

10 10 Scaling factor (M)

10

2

10

0 −1 10

3

d 15

10

10

10

10

5 10

20

30

Foxp3 total (Ft )

40

5 0 0

10

20

30

Foxp3 total (Ft )

40

Potential

c 15 Potential

b 15 Potential

a 15

0 0

H1 H2 H1*H2

Barrier Height

Foxp3 total (Ft )

24

Potential

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


8

Page 31 of 32

5 0 0

10

20

30

Foxp3 total (Ft )

40

0

1

10 10 Scaling factor (M)

10

2

5 0 0

10

20

30

Foxp3 total (Ft )

40

Figure 6: Network robustness with respect to FoxP3-DNA binding affinity. (A) Treg network transits from monostably high (Treg state) to bistable and to monostably low (Tconv state) with proportional increase of all of the dissociation constants of FoxP3-DNA complexes. The four potential landscapes (a, b, c, and d) represent different dynamic characteristics. (B) Network √ robustness ( H1 H2 ) is altered upon the scaling of the dissociation constants.



45

8

+

Robustness metric

Up-regulated TF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

30

15

Page 32 of 32

Treg

0 0

Tconv 20 40 Down-regulated TF

60

H1 H2 H1*H2

9 6 3 0

A

B

C

D


E

A Minimal Transcriptional Controlling Network of Regulatory T Cell

Recommend Documents