Modeling the dynamics of transcriptional gene regulatory networks for animal development

Abstract

The dynamic process of cell fate specification is regulated by networks of regulatory genes. The architecture of the network defines the temporal order of specification events. To understand the dynamic control of the developmental process, the kinetics of mRNA and protein synthesis and the response of the cis-regulatory modules to transcription factor concentration must be considered. Here we review mathematical models for mRNA and protein synthesis kinetics which are based on experimental measurements of the rates of the relevant processes. The model comprises the response functions of cis-regulatory modules to their transcription factor inputs, by incorporating binding site occupancy and its dependence on biologically measurable quantities. We use this model to simulate gene expression, to distinguish between cis-regulatory execution of "AND" and "OR" logic functions, rationalize the oscillatory behavior of certain transcriptional auto-repressors and to show how linked subcircuits can be dealt with. Model simulations display the effects of mutation of binding sites, or perturbation of upstream gene expression. The model is a generally useful tool for understanding gene regulation and the dynamics of cell fate specification.

Fig. 1

Processes involved in transcription and translation. (A) For mRNA synthesis these are transcription initiation, RNA polymerase translocation, mRNA processing and mRNA export from the nucleus to the cytoplasm. The level of mRNA depends also on mRNA turnover rate. The processes that control protein level are translation and protein turnover rates. (B) The initiation rate controls the number of transcripts that are generated within a given time interval. The higher is the initiation rate the more mRNA copies are produced.

Fig. 2

mRNA and protein accumulation functions. (A) The mRNA and protein accumulation curves were obtained by substituting the following kinetic parameters in eqs. (1) and (2): I_s=3 initiations/minute, k_t=2 protein/(mRNA×minute), k_dm=Ln2/60=0.012 min⁻¹, and k_dp=Ln2/40=0.014 min⁻¹. The initiation rate, I_s, is the initial linear slope of the mRNA accumulation curve. The maximal level of mRNA is the steady state level, I_s/k_dm. The half-life, t_1/2=ln2/k_dm, is the time when the mRNA accumulation function reaches half maximum. The simulation was done using Mathematica 5.2. (B) Use of the model to fit experimental measurements of mRNA levels in sea urchin embryos (Howard-Ashby et al., 2006). Left, for maternal transcripts, i.e., genes that their mRNA is present in the egg, the initiation rate is zero at early times and the mRNA level decays exponentially as e^−kdmt, Eq. (1). A fit to the measured mRNA time course for the maternal phase of the gene oct1.2, results in half-lifetime of 4.18 h. (oct1.2 has a zygotic phase, i.e., transcription that starts after fertilization, initiated at 18 hpf, that was not fitted with the model). Right, Eq. (3) was used to fit the mRNA accumulation curve for the zygotic gene, tgif. The zygotic expression of this gene starts at about 18 hpf, so 18 hpf is the t=0 in this simulation. The result is an initiation rate of 124 molecules/hour and half-life of 14.6 h. Reprinted from Howard-Ashby et al., 2006. Dev. Biol. 300, 74–89; copyright Elsevier, Inc.

Fig. 3

Occupancy and transcription. (A) The occupancy of a binding site depends on the ratio of its association and dissociation rate constants, k_aS and k_dS, respectively, and on the transcription factor concentration. (B) Initiation rate dependence on occupancy for different activation strengths, k_b (Eq. (12)). Red curve, k_b=5, green curve, k_b=20 and blue curve k_b=50. At low occupancy the initiation rate increases linearly with the occupancy with a slope of k_b. In this example we consider I_max=11 initiations per minute, as calculated in text for 2 gene copies at 15 °C. The simulation was done using Mathematica 5.2. (C) Cooperative binding to the DNA increases the stability of the factor–DNA complex. The cooperativity constant, K_q, indicates how much the two factor–DNA complex is stabilized compared to independent binding of the two factors. Free energy contributions for DNA–protein and protein–protein interactions are indicated by green and yellow arrows respectively.

Fig. 4

Simple GRN subcircuits and kinetic outputs. (A) Subcircuit in which regulatory genes a and b produce factors that activate the expression of gene c. (B) Time courses for expression of a, b and c, assuming different logic gates. Upper panel: Time courses for protein output of a (magenta) and b (cyan). b is activated 60 min after the activation of a and both factors are activated at constant initiation rate of I_s=2. Bottom panel: Time course for c mRNA under different cis-regulatory gates processing inputs from a and b genes. Red curve: c is regulated by a Additive OR b inputs, Eq. (14). Blue curve: c is regulated by a AND b inputs, Eq. (16), K_q=1. Green curve: c is regulated by a AND b inputs and the binding of a and b is cooperative, Eq. (16), K_q=20. The parameters used in this simulation are: Relative equilibrium constant, K_r=10⁵, activation strength, k_b=5, mRNA turnover rate k_dm=0.001 min⁻¹, protein turnover k_dp=0.002 min⁻¹, translation rate, k_t=2 protein/(mRNA×minute), mRNA transcription delay, T_m=20 min. The number of non-specific sites, D_n, was estimated as 90% of the total sea urchin genome, which is 8×10⁸, so D_n=7.2×10⁸. The initial levels of all genes, a, b and c was assumed to be zero at time zero. (C) Auto-repression sub-circuit. (D) Time courses of mRNA (left) and protein (right) for an auto-repressor operating according to the threshold model (Eq. 18). The kinetic parameters used in this simulation are: k_t=2, I_s=2, K_r=10⁵, k_dm=k_dp=0.017 min⁻¹, Y₀=0.36, B₀=0.2, D_n=7.2×10⁸ and T_m=20 min. Simulations were done using Mathematica 5.2.

Fig. 5

Compound GRN circuit. (A) Schematic diagram of the circuit. Gene A activates gene B. Gene B has a positive feedback into its own cis-regulatory module. Gene B activates gene C, and genes B and C together activate gene D. The cis-regulatory module of B executes Additive OR logic on A and B, and the cis-regulatory module of D executes AND logic on B and C. (B) Time courses of the mRNA expression levels of genes A (magenta), B (cyan), C (dark blue) and D (green). (C) Time course of the protein expression levels of genes A, B, C and D, color code similar to (B). The parameters used in this simulation are: Relative equilibrium constant, K_r=10⁵, activation strength for all the equations, k_b=2, mRNA turnover rates: k_dmA=0.001 min⁻¹, k_dmB=k_dmC=k_dmD= 0.005 min⁻¹, protein turnover rates: k_dp=0.01 min⁻¹, k_dpB=k_dpC=k_dpD=0.008 min⁻¹, translation rate, k_t=2 protein/(mRNA×minute), mRNA transcription delay, T_m=40 min, cooperativity factor K_q=1. The number of non-specific sites, D_n, was estimated as 90% of the total sea urchin genome, which is 8×10⁸, so D_n=7.2×10⁸. The initial levels of the protein A, and the mRNA and protein of B, C and D were assumed to be zero at time zero. The initial mRNA level of gene A, mA(0)=500 molecules per cell. Simulations were done using Mathematica 5.2.