Sensitivity Analysis of Retrovirus HTLV-1 Transactivation

Abstract

Human T-cell leukemia virus type 1 is a human retrovirus endemic in many areas of the world. Although many studies indicated a key role of the viral protein Tax in the control of viral transcription, the mechanisms controlling HTLV-1 expression and its persistence in vivo are still poorly understood. To assess Tax effects on viral kinetics, we developed a HTLV-1 model. Two parameters that capture both its deterministic and stochastic behavior were quantified: Tax signal-to-noise ratio (SNR), which measures the effect of stochastic phenomena on Tax expression as the ratio between the protein steady-state level and the variance of the noise causing fluctuations around this value; t_1/2, a parameter representative of the duration of Tax transient expression pulses, that is, of Tax bursts due to stochastic phenomena. Sensitivity analysis indicates that the major determinant of Tax SNR is the transactivation constant, the system parameter weighting the enhancement of retrovirus transcription due to transactivation. In contrast, t_1/2 is strongly influenced by the degradation rate of the mRNA. In addition to shedding light into the mechanism of Tax transactivation, the obtained results are of potential interest for novel drug development strategies since the two parameters most affecting Tax transactivation can be experimentally tuned, e.g. by perturbing protein phosphorylation and by RNA interference.

1. Introduction

Human T-cell leukemia virus type 1 (HTLV-1) is a human retrovirus first isolated in the late 1970s. It is endemic in many areas of the world, for example, in southern Japan, in the Caribbean basin, in central Africa, in central and south America, in the Melanesian Islands in the Pacific ocean, and in the aboriginal population in Australia; 15–20 million people are estimated to be HTLV-1-infected worldwide (Lairmore and Franchini, 2007; Proietti et al., 2005). HTLV-1 can cause either a neurodegenerative disease termed Tropical Spastic Paraparesis, or Adult T-cell Leukemia/Lymphoma (Lairmore and Franchini, 2007), an aggressive neoplasm that is refractory to current therapy regimens.

To gain insight into the mechanisms controlling HTLV-1 expression, we recently developed a mathematical model of HTLV-1 (Corradin et al., 2010). In this model, the noise affecting viral gene transcription and translation was studied using Gillespie's exact stochastic simulations (Gillespie, 1976, 1977, 1992, 2007). Results indicated that stochastic phenomena can sporadically induce remarkable fluctuations of the expression of protein Tax, which is the main signal related to retrovirus HTLV-1 activation (Brauweiler et al., 1995), over the expected level corresponding to its deterministic steady state. This concept is consistent with other experimental models of different viruses where bursts in protein expression were described (Kierzek et al., 2001; Cai et al., 2006; Yu et al., 2006). Since Tax transactivates the HTLV-1 gene circuit, i.e., enhances retrovirus transcription, transient pulses in its expression suggest mechanisms of retroviral activation similar to those proposed for HIV by Weinberger et al. (2008), who experimentally observed that HIV activation is influenced by the duration of the stochastic bursts of expression of the Tat protein. In analogy to HTLV-1, Tat transactivates the HIV gene circuit inducing a positive feedback mechanism whose strength is determined by the duration of Tat bursts. The role of positive autoregulation on the HIV viral genes kinetics and on the development of pathological states was recently pointed out also by other authors (Burnett et al., 2009; Althaus and De Boer, 2010). In more general terms, the effects of stochastic phenomena on the expression of proteins were studied in depth (Thattai and van Oudenaarden, 2001; Paulsson, 2004; Kaern et al., 2005; Rao et al., 2002; Ozbudak et al., 2002), with particular attention to gene circuits characterized by positive feedback mechanisms that can induce bimodalities, which are potentially related to the existence of two distinct states, e.g. latency and activation (Friedman et al., 2006; To and Maheshri, 2010).

In the present study, we aimed at gaining a better insight into the mechanisms controlling HTLV-1 activation. Since modulating transactivation might be an effective means to control HTLV-1 activation, viral spread, and possibly, disease in infected patients, we focused on the role of the transactivator protein Tax in the absence of Rex, a key HTLV-1 regulatory protein that controls viral gene expression at the post-transcriptional level (Lairmore and Franchini, 2007); this way the effects of transactivation on Tax expression were isolated. An additional advantage of this reduced circuit is that Tax repression, performed by Rex, is avoided, thus facilitating future experimental validations of the results of our work (Lin et al., 1998). Once derived a mathematical model of the reduced gene circuit from the recently developed full-length HTLV-1 model (Corradin et al., 2010), an analytic approach based on transfer functions was pursued to assess the propagation of noise affecting genes/proteins synthesis and degradation to Tax expression by means of two parameters: Tax signal-to-noise ratio (SNR), which is the ratio between the protein steady-state level and the variance of the noise causing fluctuations around this value and t_1/2, derived from the noise autocorrelation function (Priestley, 1982) and related to the duration of Tax bursts (Weinberger, 2008).

To test the reliability of the transfer function approach, originally proposed to evaluate the effects of stochastic phenomena on simple gene circuits at steady state (Simpson et al., 2003; Cox et al., 2006), results were compared to those obtained using Gillespie algorithm, which represents an acceptable in silico surrogate of biological validation.

Finally, sensitivity analysis (Thomaseth and Cobelli, 1999) was applied to identify the model parameters most affecting either SNR or t_1/2 and, thus, Tax expression and transactivation.

2. Methods

2.1. Model of the HTLV-1-derived gene circuit

The main biological processes characterizing the HTLV-1 gene circuit, obtained by depriving of Rex synthesis and effects a recently developed full-length model (Corradin et al., 2010), are schematized in Figure 1 as:

FIG. 1.

HTLV-1-derived gene circuit. N₁, N₂, N′₃, N″₃ represent the noise affecting gag synthesis, tax/rex degradation, Tax synthesis, and degradation, respectively. Parameters k_s, k₀₂, and k₀₃ represent the rate constant of gag splicing, tax/rex, and Tax degradation, respectively; β′₃ is the product of Tax gain in translation multiplied for the rate constant of the chemical reaction representing Tax translation. Continuous lines indicate fluxes; dashed lines indicate controls.

i. Synthesis of the primary transcript gag;

ii. Gag splicing into transcript tax/rex and tax/rex degradation;

iii. Translation and degradation of protein Tax, which is encoded by tax/rex and enhances gag transcription (transactivation).

The kinetics of the primary transcript gag, gene tax/rex, and protein Tax can be described by the following set of differential equations (Eqs. 1–3), which are based on mass action kinetics (Corradin et al., 2010; Connors, 1991; Kaern et al., 2003): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{cases} \frac {{\rm dq}_1} {{\rm dt}} = {\rm S}_1 ({\rm t}) - {\rm k}_{\rm s} \cdot {\rm q}_1 ({\rm t}) = {\rm MOI} \cdot {\rm S} \cdot \left\{1 + {\rm \beta}^{\prime}_1 \cdot \frac {{\rm q}_3 ({\rm t})^2} {{\rm h}^{\prime 2}_1 + {\rm q}_2 ({\rm t})^2} \right\} - {\rm k}_{\rm s} \cdot {\rm q}_1 ({\rm t}) \qquad {\rm q}_1 (0) = {\rm q}_1^{\rm SS} & (1) \\ \frac {{\rm dq}_2} {{\rm dt}} = {\rm S}_2 ({\rm t}) - {\rm k}_{02} \cdot {\rm q}_2 ({\rm t}) = {\rm k}_{\rm s} \cdot {\rm q}_1 ({\rm t}) - {\rm k}_{02} \cdot {\rm q}_2 ({\rm t}) \qquad {\rm q}_2 (0) = {\rm q}_2^{\rm SS} & (2) \\ \frac {{\rm dq}_3} {{\rm dt}} = {\rm S}_3 ({\rm t}) - {\rm k}_{03} \cdot {\rm q}_3 ({\rm t}) = {\rm \beta}_3^{\prime} \cdot {\rm q}_2 ({\rm t}) - {\rm k}_{03} \cdot {\rm q}_3 ({\rm t}) \qquad {\rm q}_3 (0) = {\rm q}_3^{\rm SS}& (3) \end{cases} \end{align*} \end{document}

where the state variables q_i, i = 1 … 3, represent the number of molecules in a cell for gag, tax/rex and Tax, respectively, and superscript SS indicates their steady-state values as initial conditions. S_i, i = 1 … 3, are synthesis processes: transcript gag synthesis, S₁, is due to a minimal transcription, MOI·S, and is enhanced by the transactivation mechanism due to protein Tax, as expressed by a nonlinear Hill term; tax/rex synthesis, S₂, corresponds to the amount of spliced gag, k_s· q₁; Tax synthesis, S₃, is proportional to tax/rex through β′₃. Parameter k_s [1/h] represents the nuclear export rate whereas k₀₂ and k₀₃ [1/h] are the degradation rates of transcript tax/rex and protein Tax, respectively. Parameter MOI represents the multiplicity of infection, i.e., the number of viral genomes integrated in the host cell genome while S [1/h] is the minimal transcription rate. As concerns the Hill term in Eq. (1), β′₁ is dimensionless and called transactivation constant, while the Hill constant h’₁ is the product of many chemical equilibrium constants, as detailed in Corradin et al. (2010). Finally, parameter β′₃ [protein molecules/(transcript molecules*h)] is the product of Tax gain in translation multiplied for the rate constant of the chemical reaction representing Tax translation (Corradin et al., 2010).

Since the basal transcription rate S and Hill constant h′₁ are unknown, model equations were scaled for these two parameters, leading to the substitution of time t with the dimensionless τ = t · S (Hasty et al., 2000, 2001), and to the normalization of the state variable q₃, which was substituted by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\widetilde {\rm q}_3} = {\rm q}_3 / {{\rm h}^{\prime}}_{1}$$\end{document} . As a consequence, model equations became: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{cases} \frac {{\rm dq}_1} {{\rm d \tau}} = {\rm MOI} \cdot \left\{1 + {\rm \beta}^{\prime}_1 \cdot \frac {{\widetilde {\rm q}_3} ({\rm \tau})^2} {1 + {\widetilde {\rm q}_3} ({\rm \tau})^2} \right\} - {\rm k}_{\rm s}^{\prime} \cdot {\rm q}_1 ({\rm \tau}) \qquad {\rm q}_1 (0) = {\rm q}_1^{\rm SS} & (4) \\ \frac {{\rm dq}_2} {{\rm d \tau}} = {\rm k}_{\rm s}^{\prime} \cdot {\rm q}_1 ({\rm \tau}) - {\rm k}_{02}^{\prime} \cdot {\rm q}_2 ({\rm \tau}) \qquad {\rm q}_2 (0) = {\rm q}_2^{\rm SS} & (5) \\ \frac {{\rm d} \widetilde {\rm q}_3} {{\rm d \tau}} = \widetilde {\rm \beta}_3^{\prime} \cdot {\rm q}_2 ({\rm \tau}) - {\rm k}^{\prime}_{03} \cdot \widetilde {\rm q}_3 ({\rm \tau}) \qquad \widetilde {\rm q}_3 (0) = \widetilde {\rm q}_3^{\rm SS} & (6) \end{cases} \end{align*} \end{document}

Nominal values for all the parameters were fixed according to Brauweiler et al. (1995): β′₁ = 10, k′s = k_s/S = 1, k′₀₂ = k₀₂/S = 0.01, k′₀₃ = k₀₃/S = 0.01, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm \beta}^\prime_3 = 0.01$$\end{document} (dimensionless), and MOI = 1 viral particle/cell.

2.2. Modeling intrinsic noise

The biological processes of transcript and protein synthesis and degradation are affected by intrinsic noise (Thattai and van Oudenaarden, 2001; Kaern et al., 2005). Thus, noise sources N₁, N₂, N′_3, and N′′₃, were added to the model (Fig. 1) to represent the noise affecting gag transcription, tax/rex degradation, Tax synthesis, and degradation, respectively. As in Simpson et al. (2003), the stochastic phenomena affecting transcription and translation were modeled as shot noise, which characterizes many electronic and optical devices when the number of particles that carry energy (electrons or photons) is finite and small. This description appears to be reasonable also for biological noise, in particular for the noise affecting transcription since the number of retrovirus promoters in HTLV-1-infected patients is usually very limited (Taylor et al., 1999).

From a mathematical point of view, shot noise is a Poisson process where events occur independently; therefore, it is characterized, in the frequency domain by a uniform distribution of power (Priestley, 1982), i.e., its power spectral density (PSD) is constant for all positive frequencies and equal to either molecular synthesis (as regards the noise in transcription and translation) or molecular degradation (as regards the noise in degradation) multiplied by two (Simpson et al., 2003). Therefore, PSDs of noise sources N₁ (affecting gag synthesis) and N₂ (affecting tax/rex degradation), named PSD₁ and PSD₂, respectively, do not vary with frequency f and can be expressed as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm PSD}_1 ({\rm f}) & = 2 \cdot {\rm S}_1^{\rm SS} = 2 \cdot {\rm MOI} \cdot \left\{1 + {\rm \beta}^{\prime}_1 \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} \right\} & (7) \\ {\rm PSD}_2 ({\rm f}) & = 2 \cdot {\rm S}_2^{\rm SS} = 2 \cdot {{\rm k}^{\prime}}_{02} \cdot {{\rm q}_2}^{\rm SS} = 2 \cdot {\rm MOI} \cdot \left\{1 + {\rm \beta}^{\prime}_1 \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} \right\} & (8) \end{align*} \end{document}

where S₁^SS and S₂^SS represent the syntheses of gag and tax/rex at steady state, respectively, obtained by solving Eqs. (4–6) at steady state (see Appendix). Following the definition of shot noise, since at equilibrium the degradation of protein Tax equals its synthesis, noise terms N′₃ and N′′₃ have equal PSDs. Assuming they are independent, they can be considered as a unique noise source, denoted as N₃, with constant PSD, PSD₃, equal to the sum of the PSDs of N′₃ and N′′₃ (Simpson et al., 2003), and given by: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm PSD}_3 ({\rm f}) = 2 \cdot 2 \cdot {\rm S}_3^{\rm SS} = 2 \cdot 2 \cdot {\widetilde {\rm \beta}^{\prime}}_{3} \cdot {{\rm q}_2}^{\rm SS} = 2 \cdot 2 \cdot {\widetilde {\rm \beta}^{\prime}}_{3} \cdot \frac {{\rm MOI}} {{\rm k}^{\prime}_{02}} \cdot \left\{ 1 + {\rm \beta}^{\prime}_1 \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} \right\} \tag {9} \end{align*} \end{document}

where S₃^SS represents Tax synthesis at steady state.

Stochastic phenomena affecting synthesis and degradation of gag, tax/rex, and Tax induce random deviations of Tax from its steady-state value, modeled as an additive noise N_out (Fig. 1). Its power spectral density, PSD_out, can be calculated from PSD₁, PSD₂ and PSD₃ in Eqs. (7–9) by means of the transfer function approach: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm PSD}_{\rm out} = {\rm PSD}_1 ({\rm f}) \cdot \mid {\rm H}_{41} (\rm f) \mid^2 + {\rm PSD}_2 ({\rm f}) \cdot \mid {\rm H}_{42} ({\rm f}) \mid^2 + {\rm PSD}_3 ({\rm f}) \cdot \mid {\rm H}_{43} ({\rm f}) \mid^2 \tag{10} \end{align*} \end{document}

where f is the frequency and H_4i, I = 1 … 3 are related to the transfer functions of the system in Eqs. (4–6), expressing the relationship between N_out and the noise sources N_i, I = 1 … 3 as a function of model parameters, as detailed in the Appendix. Since this approach is applicable only to linear systems while the model in Eqs. (4–6) is nonlinear (because of the Hill term in Eq. 4), a linearization around a steady-state solution was performed, as detailed in the Appendix. By applying the inverse Fourier transform to PSD_out, the autocorrelation function ACF_out can be analytically calculated (Priestley, 1982) (details in the Appendix), which provides important information on the relationship among noise samples taken at different. In fact, ACF_out(k) represents the correlation among samples separated by the time interval k and its value for k = 0 represents the standard deviation of the noise N_out, denoted as σ_out (Priestley, 1982), from which SNR can be calculated: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm SNR} = \frac {\widetilde {\rm q}_3^{\rm SS}} {\sigma_{\rm out}} = \frac {\widetilde {\rm q}_3^{\rm SS}} {\sqrt {{\rm ACF}_{\rm out} (0)}} \tag {11} \end{align*} \end{document}

Moreover, as suggested by Weinberger et al. (2008), the time lag corresponding to the halving of ACF_out, denoted as t_1/2, is indicative of transient pulses duration, and can be a useful parameter since longer bursts in the expression of the transactivator protein are likely to favor retrovirus activation (Weinberger et al., 2008).

2.3. Sensitivity analysis

In order to identify the model parameters most affecting the gene circuit SNR, sensitivity analysis was applied (Thomaseth and Cobelli, 1999), i.e., sensitivity of SNR to model parameter p_i, S_i, was defined as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm S}_{\rm i} = \frac {\rm {dSNR} ({\rm p}_1 , \ldots , {\rm p}_{\rm n})} {{\rm dp}_{\rm i}} \cdot \frac {{\rm p}_{\rm i}} {\rm {SNR} ({\rm p}_1 , \ldots , {\rm p}_{\rm n})} \tag {12} \end{align*} \end{document}

and evaluated from the analytical expression of the SNR provided by Eq. (11).

Since it was not possible to derive an analytical expression of t_1/2 from that of ACF_out, we investigated the effects exerted on t_1/2 by model parameters in terms of Spearman correlation. In more detail, the analytic expression of ACF_out was evaluated for different values of one model parameter at a time: for each ACF_out time course, parameter t_1/2 was estimated and Spearman correlation between the set of model parameter values and the set of extracted t_1/2 values was calculated.

2.4. Gillespie's exact stochastic simulations

Tax steady-state fluctuations were simulated by Dizzy software package (Ramsey et al., 2005), which implements the direct method of Gillespie algorithm (Gillespie, 1976, 1977, 1992, 2007). 100 distinct simulations were performed; from each of them noise samples were extracted following the de-trending procedure described in Weinberger et al. (2008). A mean noise autocorrelation function was estimated from the whole set of 100 simulations using the biased estimator: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm ACF}_{\rm out} ({\rm k}) = \frac {\sum \nolimits_{{\rm m} = 1}^{\rm M} \sum \nolimits_{{\rm n} = 1}^{{\rm N} - {\rm k} - 1} \widetilde {\rm N}_{\rm m} ({\rm n} + {\rm k}) \widetilde {\rm N}_{\rm m} ({\rm n})} {{\rm N \cdot M}} \tag {13} \end{align*} \end{document}

where k = 0 … N−1 is the scaled time lag, subscript m = 1 … 100 refers to the m^th simulation; thus, N_m(n) are the noise samples extracted from the m^th simulation and N = 500,001 is the number of time samples composing each time series. Hence, numerical PSD_out was estimated via discrete Fourier transform algorithm (Smith, 2007).

3. Results

The time courses of PSD_out (Eq. 10) and ACF_out are shown, plotted in continuous line, in Figure 2 where the x-axes are scaled frequency or time lag because of the time scaling performed on the model. PSD and ACF numerically estimated from Gillespie's steady-state fluctuations (Gillespie, 1976, 1977, 1992, 2007) are superimposed, plotted in dashed line. Their close match confirms the reliability of the noise propagation procedure based on transfer functions to analyze HTLV-1 gene circuit. Both PSDs are characterized by a relevant fall, which indicates that the noise affecting Tax expression is characterized mainly by low frequencies, at variance with the original noise sources N₁, N₂ and N₃, assumed to have constant PSDs and thus uniform frequency content. This fall is due to the presence of zeros at the denominator (poles) of the transfer functions H_4i, i = 1 … 3 and is expected since it is well known that the presence of a pole in the denominator of a transfer function causes a decrease of 20 dB/decade in its magnitude (Levine, 1996). In more detail, the fall in the trend of PSD_out is due to the combined effect of two close poles in the denominators of H_4i, which are located at the scaled frequencies 1.5e-3 and 1.7e-3, respectively. These values correspond to 2.9e-6 Hz and 3.1e-6 Hz, if the scaling factor S is set to adenovirus value (Kugel and Goodrich, 1998), which is 1.9e-3/s.

FIG. 2.

PSD and ACF of the noise affecting Tax expression. Analytic PSD_out (a; solid black line) and ACF_out (b; solid black line). Their counterparts, estimated from Gillespie's exact stochastic simulations, are superimposed in dashed line. Values of parameter t_1/2, resulting from the two approaches, are indicated along the x-axis.

Numerical values of the gene circuit SNR, of Tax steady-state expression level and of the variance characterizing noise N_out are shown in Table 1, together with the time lag t_1/2 at which ACF_out assumes a value equal to 50% of the maximum level. The corresponding quantities, numerically estimated from Gillespie's steady-state fluctuations, are reported in Table 1, showing a good agreement between the two methods. In particular, t_1/2 values, estimated via either the transfer function approach or Gillespie's exact stochastic simulations, equal 99 and 115, respectively, when expressed in scaled time, which correspond to 14.5 and 16.8 h if the scaling factor S is set as before.

Table 1.

Analytic Results Compared to Numerical Estimates

	Transfer function approach	Gillespie's exact stochastic simulations
SNR	24.88	23.66
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm q}^{\rm SS}_3$$\end{document}	1,100 molecules	1,100 molecules
N_out variance	1,954.34 molecules²	2,162.20 molecules²
t_1/2	99	115

Values of Tax SNR and steady-state expression \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm q}^{\rm SS}_3$$\end{document} , variance of the noise N_out affecting Tax, and the time lag t_1/2 corresponding to the halving of ACF_out, estimated via either the transfer function approach or Gillespie's exact stochastic simulations.

Results of sensitivity analysis, summarized in Table 2, suggest a major role of MOI in increasing SNR and an opposite role of Tax degradation rate, k′₀₃. Parameter \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm \beta}^{\prime}_1$$\end{document} , which represents the transactivation constant, is also a relevant determinant of SNR while parameter \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm \beta}^{\prime}_3$$\end{document} , related to Tax translation, plays a minor role in increasing SNR, equal and opposite to that of tax/rex degradation rate, k′₀₂. Parameter k′_s does not affect SNR, since it has no role neither on synthesis nor on degradation, which are the biological processes determining Tax presence and permanence in the cell and, consequently, its steady-state level and SNR.

Table 2.

SNR Sensitivities and t_1/2 Spearman Correlations

	MOI	β′₁	k′₀₂	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm \beta}^{\prime}_3$$\end{document}	k′₀₃	k′_s
SNR sensitivities	0.500	0.454	−0.281	0.281	−0.500	0
t_1/2 correlations	0	0	−0.946 (p < 0.01)	1 (p < 0.01)	−0.946 (p < 0.01)	0

Influence of model parameters on SNR and t_1/2, measured from sensitivity analysis for the former and Spearman correlation for the latter.

Spearman correlation coefficients suggest that degradation rates of transcript tax/rex, k′₀₂, and protein Tax, k′₀₃, negatively correlate with t_1/2, which is reasonable because an increased degradation means a shorter permanence of the protein in the cell. Equal and opposite effect is exerted by Tax gain in translation, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm \beta}^{\prime}_3$$\end{document} , able to determine an increase in Tax synthesis. On the other hand, the other system parameters, MOI, β′₁, and k′_s, do not correlate significantly with parameter t_1/2.

4. Conclusion

In this article, we focused on the biological mechanism of transactivation, which consists in the enhancement of viral transcription; this is due to the HTLV-1 protein Tax that, consequently, is thought to be related to retrovirus activation. Such a relationship was recently observed in other biological models concerning retrovirus HIV and its transactivator protein Tat.

In the present work, in order to get a better insight on HTLV-1 activation, we characterized Tax expression on a HTLV-1-derived gene circuit, appropriately modified to isolate the effects of the transactivation mechanism. The mathematical model of this HTLV-1-derived gene circuit was obtained from a recently developed full-length HTLV-1 model (Corradin et al., 2010), and includes the kinetics of transcripts gag and tax/rex, and of protein Tax. An analytic approach based on transfer functions was adopted to assess the propagation of noise affecting gene/protein synthesis and degradation to Tax expression and to estimate two parameters able to quantify both its deterministic and stochastic behavior: (i) the HTLV-1 gene circuit SNR, i.e., Tax steady-state level normalized to the standard deviation of the noise affecting this protein expression and (ii) t_1/2, a parameter representative of the duration of Tax transient expression pulses, i.e., of Tax bursts due to stochastic phenomena that can cause important deviations of the protein expression from its steady-state value.

As concerns the assumption of shot noise, its limits to model stochastic phenomena in biological systems were discussed by Simpson et al. (2003); moreover, Golding et al. (2005) recently observed mRNA distributions larger than the Poissonian, which underlies the shot noise. However, also, the Gillespie algorithm is based on Possonian statistics (Wilkinson, 2006), and Poissonian noise was recently re-proposed by Thattai and van Oudenaarden (2001). As concerns the other two assumptions that characterize shot noise, i.e., the independence of noise samples and the constant power spectral density equal to the molecular turn-over multiplied by two, they were not experimentally proved but the agreement with the results derived from Gillespie's exact stochastic simulations suggest that the hypotheses underlying shot noise are reasonable in the case of HTLV-1 model presented in this work. Similar good agreements were reported also for different gene circuits characterized by negative feedback loops (Simpson et al., 2003, 2004).

On the other hand, Gillespie's steady-state fluctuations are much more time-consuming than the approach based on shot noise and transfer functions; moreover, several computational steps are needed to extract noise samples from the time series of protein Tax expression in the de-trending procedure (Weinberger et al., 2008). Here, we obtained similar results by means of a much simpler and immediate approach that can be particularly helpful in a synthetic biology context where it is of interest to evaluate noise propagation. Moreover, the description of noise in the domain of frequency before performing a biological experiment allows introducing noise rejection techniques in the design of gene circuits.

Given such analytic framework, sensitivity analysis was performed. Some results are of potential interest in drug target selection since two model parameters most able to influence SNR and t_1/2 were identified: transactivation constant and tax/rex degradation rate. These two parameters can be experimentally tuned: transactivation constant can be affected by preventing Tax phosphorylation with appropriate site-directed mutagenesis (Bex et al., 1999) or, on the opposite, by enhancing this biological process as described in Fontes et al. (1993), and tax/rex degradation rate can be remarkably increased by RNA silencing.

Future developments concern the experimental validation of the main results reported in this article; in particular, testing (i) if affecting Tax phosphorylation, the SNR actually changes letting unaltered Tax expression bursts (Table 2) and (ii) if an increase in tax/rex degradation rate reflects in the shortening of the duration of protein Tax transient expression pulses decreasing limitedly gene circuit SNR (as indicated in Fig. 3 and Table 2). These validations should be feasible since the HTLV-1-derived gene circuit schematized in Fig. 1 can be easily implemented in laboratory by deletion of protein Rex initiation codon, as done in Gröne et al. (1996).

FIG. 3.

Relationship between protein Tax bursts and tax/rex degradation rate. Realizations of Gillespie's steady-state fluctuations for different values of tax/rex degradation rate k′₀₂ show that the amplitude of protein Tax bursts markedly reduce due to an increase of k′₀₂.

As regards the theoretical analysis of the retrovirus gene circuit, it can be extended by considering different noise models, e.g., based on the experimental observations of Golding et al. (2005). A further extension regards modeling the extrinsic noise (Elowitz et al., 2005), i.e., the variability due to different amounts of cellular molecules like RNA polymerase or cellular transcription factors. Here, this noise component was neglected because, with low multiplicities of infection, these reactants are in great excess with respect to the viral promoters they interact with; hence, they behave like constant amounts in the chemical kinetics (Connors, 1991) characterizing retrovirus transcription. However, it can be of interest to address the higher MOI case, together with modeling and quantification of extrinsic noise.

In conclusion, in this article we identified two biological processes (Tax phosphorylation and tax/rex degradation) that potentially affect remarkably transactivation, but with distinct effects on the expression of the transactivator protein Tax. This offers some hypotheses, e.g., related to drug target selection, to affect retrovirus expression and the duration of the asymptomatic state in HTLV-1-infected patients, to be experimentally validated.

5. Appendix

5.1. Molecular syntheses and degradations at steady state

Following the definition of shot noise, PSD₁ is constant and equal to gag synthesis at steady state, S₁^SS, multiplied by two (Simpson et al., 2003). S₁^SS can be obtained by solving Eqs. (4–6) at steady state: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm S}_1^{\rm SS} = {\rm MOI} \cdot \left\{1 + {\rm \beta}_1^{\prime} \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} \right\} \tag {14} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm q}^{\rm SS}_3$$\end{document} represents the steady-state expression value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm q}_3$$\end{document} . Analogously, PSD₂ is constant and equal to tax/rex synthesis, S₂^SS, multiplied by two (Simpson et al., 2003). Since synthesis and degradation are equal at steady state, S₂^SS can be calculated as k′₀₂· q₂^SS (see Eq. 5), where q₂^SS represents the steady-state expression value of q₂: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm S}_2^{\rm SS} = {\rm MOI} \cdot \left\{1 + {\rm \beta}_1^{\prime} \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} \right\} \tag {15} \end{align*} \end{document}

Similarly, synthesis and degradation of protein Tax are equal at steady state. Tax synthesis at steady state, S₃^SS, can be calculated as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm \beta}^{\prime}_3 \cdot {\rm q}_2^{\rm SS}$$\end{document} (see Eq. 6): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm S}_3^{\rm SS} = \widetilde {\rm \beta}^{\prime}_3 \cdot {\rm q}_2^{\rm SS} = \widetilde {\rm \beta}^{\prime}_3 \cdot \frac {{\rm MOI}} {{\rm k}^{\prime}_{02}} \cdot \left\{1 + {\rm \beta}_1^{\prime} \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} \right\} \tag {16} \end{align*} \end{document}

5.2. HTLV-1 gene circuit transfer functions

Since the transfer function approach can be applied only to linear systems, linearization of the HTLV-1 model (Eqs. 4–;6) around its steady state is needed. Linearization is justified a posteriori by the small coefficient of variation (4%), the inverse of the gene circuit SNR (24.88; Table 1), which confirms that the fluctuations of Tax expression due to stochastic phenomena are small, on average, with respect to the protein steady-state level.

The unique nonlinear term in the model is the Hill term of Eq. (4) whose Taylor expansion is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac {\widetilde {{\rm q}_3}^2} {1 + \widetilde {{\rm q}_3}^2} \approx \frac {\widetilde {{\rm q}_3}^{{\rm SS}^2}} {1 + \widetilde {{\rm q}_3}^{{\rm SS}^2}} + \frac {2 \widetilde {{\rm q}_3}^{{\rm SS}^2}} {[1 + {(\widetilde {{\rm q}_3}^{\rm SS})}^2]^2} \cdot (\widetilde {\rm q}_3 - \widetilde {\rm q}_3^{\rm SS}) \tag {17} \end{align*} \end{document}

where the higher order terms in the Taylor series have been neglected because we assume only small excursions around \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm q}^{\rm SS}_3$$\end{document} , which is the steady-state value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde {\rm q}_3$$\end{document} . By substituting the nonlinear Hill term with its linearized analogous in Eq. (17), model equations become: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{cases} \frac {{\rm dq}_1} {{\rm d \tau}} = {\rm MOI} \cdot {\rm \beta}^{\prime}_1 \cdot \frac {2 \widetilde {{\rm q}_3}^{{\rm SS}^2}} {\Big[1 + {\big(\widetilde {{\rm q}_3}^{\rm SS} \big)}^2 \Big]^2} \cdot \widetilde {\rm q}_3 - {\rm k}_{\rm s}^{\prime} \cdot {\rm q}_1 ({\rm \tau}) + {\rm u} \qquad {\rm q}_1 (0) = {\rm q}_1^{\rm SS} & (18) \\ \frac {{\rm dq}_2} {{\rm d \tau}} = {\rm k}_{\rm s}^{\prime} \cdot {\rm q}_1 ({\rm \tau}) - {\rm k}_{02}^{\prime} \cdot {\rm q}_2 ({\rm \tau}) \qquad {\rm q}_2 (0) = {\rm q}_2^{\rm SS} & (19) \\ \frac {{\rm d} \widetilde {\rm q}_3} {{\rm d \tau}} = \widetilde {\rm \beta}_3^{\prime} \cdot {\rm q}_2 ({\rm \tau}) - {\rm k}^{\prime}_{03} \cdot \widetilde {\rm q}_3 ({\rm \tau}) \qquad \widetilde {\rm q}_3 (0) = \widetilde {\rm q}_3^{\rm SS} & (20) \end{cases} \end{align*} \end{document}

where u represents a constant transcriptional input of the linearized model in Eqs. (18–20): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm u} = {\rm MOI} \cdot \left\{1 + {\rm \beta}_1^{\prime} \cdot \frac {(\widetilde {\rm q}_3^{\rm SS})^2} {1 + (\widetilde {\rm q}_3^{\rm SS})^2} - {\rm \beta}^{\prime}_1 \cdot \frac {2 \widetilde {{\rm q}_3}^{{\rm SS}}} {[1 + {(\widetilde {{\rm q}_3}^{\rm SS})}^2]^2} \cdot \widetilde {\rm q}_3^{\rm SS} \right\} \tag {21} \end{align*} \end{document}

A Langevin approach (Simpson et al., 2003, 2004) is then applied: additive noise terms η_i(τ), i = 1 … 3, corresponding to noise sources N_i, i = 1 … 3, (Fig. 1) and related to the noise in gag transcription, tax/rex degradation, Tax synthesis and degradation, respectively, are introduced into the linearized HTLV-1 model of Eqs. (18–20): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{cases} \frac {{\rm dq}_1} {{\rm d \tau}} = {\rm MOI} \cdot {\rm \beta}^{\prime}_1 \cdot \frac {2 \widetilde {{\rm q}_3}^{{\rm SS}^2}} {[1 + {(\widetilde {{\rm q}_3}^{\rm SS})}^2]^2} \cdot \widetilde {\rm q}_3 ({\rm \tau}) - {\rm k}_{\rm s}^{\prime} \cdot {\rm q}_1 ({\rm \tau}) + {\rm u} + {\rm \eta}_1 (\rm \tau) \qquad {\rm q}_1 (0) = {\rm q}_1^{\rm SS} & (22) \\ \frac {{\rm dq}_2} {{\rm d \tau}} = {\rm k}_{\rm s}^{\prime} \cdot {\rm q}_1 ({\rm \tau}) - {\rm k}_{02}^{\prime} \cdot {\rm q}_2 ({\rm \tau}) + {\rm \eta}_2 (\rm \tau) \qquad {\rm q}_2 (0) = {\rm q}_2^{\rm SS} & (23) \\ \frac {{\rm d} \widetilde {\rm q}_3} {{\rm d \tau}} = \widetilde {\rm \beta}_3^{\prime} \cdot {\rm q}_2 ({\rm \tau}) - {\rm k}^{\prime}_{03} \cdot \widetilde {\rm q}_3 ({\rm \tau}) + {\rm \eta}_3 (\rm \tau) \qquad \widetilde {\rm q}_3 (0) = \widetilde {\rm q}_3^{\rm SS} & (24) \end{cases} \end{align*} \end{document}

Finally, the linear differential equations in Eqs. (22–24) are Laplace transformed (Levine, 1996) as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm sQ}_1 ({\rm s}) - {\rm q}_1^{\rm SS} & = {\rm G} \cdot \widetilde{\rm Q}_3 ({\rm s}) - {\rm k}^{\prime}_{\rm s} \cdot {\rm Q}_1 ({\rm s}) + \rm{U(s)} + {\rm \eta}_1 ({\rm s}) & (25) \\{\rm sQ}_2 ({\rm s}) - {\rm q}_2^{\rm SS} & = {\rm k}^{\prime}_{\rm s} \cdot {\rm Q}_1 ({\rm s}) - {\rm k}^{\prime}_{02} \cdot {\rm Q}_2 ({\rm s}) + {\rm \eta}_2 ({\rm s}) & (26) \\ {\rm s} \widetilde{\rm Q}_3 ({\rm s}) - \widetilde{\rm q}_3^{\rm SS} & = \widetilde{\rm \beta}^{\prime}_3 \ {\rm Q}_2 ({\rm s}) - {\rm k}^{\prime}_{03} \cdot \widetilde{\rm Q}_3 ({\rm s}) + {\rm \eta}_3 ({\rm s}) & (27) \end{align*} \end{document}

where Q_i, i = 1 … 3, are the Laplace transforms of the state variables q_i; U(s) is the Laplace transform of the constant transcriptional input u and η_i(s), i = 1 … 3, those of the noise terms η_i(τ). G is the gain of the feedback loop, as results from the linearization of the Hill term in Eq. (4): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm G} = {\rm MOI} \cdot {\rm \beta}^{\prime}_1 \cdot \frac {2 \widetilde {{\rm q}_3}^{{\rm SS}}} {[1 + {(\widetilde {{\rm q}_3}^{\rm SS})}^2]^2} \tag {28} \end{align*} \end{document}

The transfer functions expressing the relationships between each noise source and protein Tax expression were found by neglecting the initial values of state variables (as usual when analyzing the forced response of any linear system [Levine, 1996]), considering one noise term at a time as the unique model input and calculating the gene circuit output, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde{\rm Q}_3 ({\rm s})$$\end{document} in our case, in function of the selected input. For sake of clarity, we report below the basic algebraic passages describing how H₄₁ was calculated: Q₁(s) in Eq. (25) is expressed as a function of η₁(s), as reported in Eq. (29); this allowed expressing Q₂(s) (Eq. 26) as a function of η₁(s) (Eq. 30): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{cases} {\rm Q}_1 ({\rm s}) = {\frac {{\rm G} \cdot \widetilde {\rm Q}_3 ({\rm s})} {{\rm s}} + {\rm k}^{\prime}_{\rm s}} + \frac {{\rm \eta}_1 ({\rm s})} {{\rm s} + {\rm k}^{\prime}_{\rm s}} &(29) \\ {\rm Q}_2 ({\rm s}) = \frac {{\rm k}^{\prime}_{\rm s}} {{\rm s} + {\rm k}^{\prime}_{02}} \cdot \left\{\frac {{\rm G} \cdot \widetilde {\rm Q}_3 ({\rm s})} {{\rm s} + {\rm k}^{\prime}_{\rm s}} + \frac {{\rm \eta}_1 ({\rm s})} {{\rm s} + {\rm k}^{\prime}_{\rm s}} \right\} & (30) \\ \widetilde {\rm Q}_3 ({\rm s}) = \frac {\widetilde {\rm \beta}^{\prime}_2} {{\rm s} + {\rm k}^{\prime}_{02}} \cdot \frac {{\rm k}^{\prime}_{\rm s}} {{\rm s} + {\rm k}^{\prime}_{02}} \cdot \left\{\frac {{\rm G} \cdot \widetilde {\rm Q}_3 ({\rm s})} {{\rm s} + {\rm k}^{\prime}_{\rm s}} + \frac {{\rm \eta}_1 ({\rm s})} {{\rm s} + {\rm k}^{\prime}_{\rm s}} \right\} &(31) \end{cases} \end{align*} \end{document}

and, subsequently, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\tilde {\rm Q}_3({\rm s})$$\end{document} (Eq. 31) as a function of η₁(s): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \widetilde {\rm Q}_3 ({\rm s}) = {\rm \eta}_1 ({\rm s}) \cdot \frac {\widetilde {\rm \beta}^{\prime}_2 \cdot {\rm k}^{\prime}_{\rm s}} {({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{\rm s}) - \widetilde {\rm \beta}^{\prime}_2 \cdot {\rm k}^{\prime}_{\rm s} \cdot {\rm G}} \tag {32} \end{align*} \end{document}

Therefore, transfer function H₄₁ can be easily calculated from Eq. (32) as the ratio between output \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\tilde {\rm Q}_3({\rm s})$$\end{document} and the input η₁(s): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm H}_{41} ({\rm s}) = \frac {\widetilde {\rm Q}_3 ({\rm s})} {{\rm \eta}_1 ({\rm s})} = \frac {\widetilde {\rm \beta}^{\prime}_2 \cdot {\rm k}^{\prime}_{\rm s}} {({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{\rm s}) - \widetilde {\rm \beta}^{\prime}_2 \cdot {\rm k}^{\prime}_{\rm s} \cdot {\rm G}} \tag {33} \end{align*} \end{document}

Analogously, H₄₂ and H₄₃ were calculated as the ratio between \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\widetilde{\rm Q}_3 ({\rm s})$$\end{document} and η₂(s) and η₃(s), respectively: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm H}_{42} ({\rm s}) & = \frac {\widetilde {\rm Q}_3 ({\rm s})} {{\rm \eta}_2 ({\rm s})} = \frac {\widetilde {\rm \beta}^{\prime}_2 \cdot ({\rm s} + {\rm k}^{\prime}_{\rm s})} {({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{\rm s}) - \widetilde {\rm \beta}^{\prime}_2 \cdot {\rm k}^{\prime}_{\rm s} \cdot {\rm G}} & (34) \\ {\rm H}_{43} ({\rm s}) & = \frac {\widetilde {\rm Q}_3 ({\rm s})} {{\rm \eta}_3 ({\rm s})} = \frac {({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{\rm s})} {({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{02}) \cdot ({\rm s} + {\rm k}^{\prime}_{\rm s}) - \widetilde {\rm \beta}^{\prime}_2 \cdot {\rm k}^{\prime}_{\rm s} \cdot {\rm G}} & (35) \end{align*} \end{document}

Finally, the system transfer functions can be expressed as a function of frequency: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm H}_{41} ({\rm f}) = \mid {\rm H}_{41} ({\rm s}) \mid_{{\rm s} = 2 {\rm j \pi f}} \tag{36} \end{align*} \end{document}

where j is the imaginary unit and f the frequency (Priestley, 1982).

5.3. Analytic calculation of ACF_out

ACF_out was calculated by symbolic inverse Fourier transform of PSD_out (Eq. 10) using MatLab (www.mathworks.it/) built-in functions.

Footnotes

Acknowledgments

A.C. thanks Sara Valpione for suggestions. This work was supported by CARIPARO 2008/2010 ("Systems Biology Approaches to Infer Gene Regulation from Gene and Protein Time Series Data").

Disclosure Statement

No competing financial interests exist.

References

Althaus

C.L.

, De Boer

R.J.

2010. Intracellular transactivation of HIV can account for the decelerating decay of virus load during drug therapy. Mol. Syst. Biol., 6:348.

Bex

, Murphy

, Wattiez

et al. 1999. Phosphorylation of the human T-cell leukemia virus type 1 transactivator Tax on adjacent serine residues is critical for Tax activation. J. Virol., 73:738–745.

Brauweiler

, Garl

, Franklin

A.A.

et al. 1995. A molecular mechanism for human T-cell leukemia virus latency and Tax transactivation. J. Biol. Chem., 270:12814–12822.

Burnett

J.C.

, Miller-Jensen

, Shah

P.S.

et al. 2009. Control of stochastic gene expression by host factors at the HIV promoter. PLoS Pathog., 5:e1000260.

Cai

, Friedman

, Xie

X.S.

2006. Stochastic protein expression in individual cells at the single molecule level. Nature, 440:358–62.

Connors

K.A.

1991. Chemical Kinetics: The Study of Reaction Rates in Solution. VCH Publishers: Weinheim, Germany.

Corradin

, Di Camillo

, Rende

et al. 2010. Retrovirus HTLV-1 gene circuit: a potential oscillator for eukaryotes. Pac. Symp. Biocomput., 2010:421–432.

Cox

C.D.

, McCollum

J.M.

, Austin

D.W.

et al. 2006. Frequency domain analysis of noise in simple gene circuits. Chaos, 16:026102.

Elowitz

M.B.

, Levine

A.J.

, Siggia

E.D.

et al. 2005. Stochastic gene expression in a single cell. Science, 297:1183–1186.

10.

Fontes

J.D.

, Strawhecker

J.M.

, Bills

N.D.

et al. 1993. Phorbol esters modulate the phosphorylation of human T-cell leukemia virus type I Tax. J. Virol., 67:4436–4441.

11.

Friedman

, Cai

, Xie

X.S.

2006. Linking stochastic dynamics to population distribution: an analytical framework of gene expression. Phys. Rev. Lett., 97:168302.

12.

Gillespie

D.T.

1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22:403–434.

13.

Gillespie

D.T.

1977. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem., 81:2340–2361.

14.

Gillespie

D.T.

1992. A rigorous derivation of the chemical master equation. Phys. A Stat. Theor. Phys, 188:404–425.

15.

Gillespie

D.T.

2007. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem., 58:35–55.

16.

Golding

, Paulsson

, Zawilski

S.M.

et al. 2005. Real-time kinetics of gene activity in individual bacteria. Cell, 123:1025–1036.

17.

Gröne

, Koch

, Grassmann

1996. The HTLV-1 Rex protein induces nuclear accumulation of unspliced viral RNA by avoiding intron excision and degradation. Virology, 218:316–325.

18.

Hasty

, Isaacs

, Dolnik

et al. 2001. Designer gene networks: towards fundamental cellular control. Chaos, 11:207–220.

19.

Hasty

, Pradines

, Dolnik

et al. 2000. Noise-based switches and amplifiers for gene expression. Proc. Natl. Acad. Sci. USA, 97:2075–2080.

20.

Kaern

, Blake

W.J.

, Collins

J.J.

2003. The engineering of gene regulatory networks. Annu. Rev. Biomed. Eng, 5:179–206.

21.

Kaern

, Elston

T.C.

, Blake

W.J.

et al. 2005. Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet., 6:451–464.

22.

Kierzek

A.M.

, Zaim

, Zielenkiewicz

2001. The effect of transcription and translation initiation frequencies on the stochastic fluctuations in prokaryotic gene expression. J. Biol. Chem, 276:8165–8172.

23.

Kugel

J.F.

, Goodrich

J.A.

1998. Promoter escape limits the rate of RNA polymerase II transcription and is enhanced by TFIIE, TFIIH, and ATP on negatively supercoiled DNA. Proc. Natl. Acad. Sci. USA, 95:9232–9237.

24.

Lairmore

M.D.

, Franchini

2007. Human T-cell leukemia virus types 1 and 2, 2071–;2106. Knipe

D.M.

, Howley

P.M.

Fields Virology, 5^th. Lippincott-Raven: Philadelphia.

25.

Levine

W.S.

1996. The Control Handbook. CRC Press: Boca Raton, FL.

26.

Lin

H.C.

, Dezzutti

C.S.

, Lal

R.B.

et al. 1998. Activation of human T-cell leukemia virus type 1 Tax gene expression in chronically infected T cells. J. Virol., 72:6264–6270.

27.

Ozbudak

E.M.

, Thattai

, Kurtser

et al. 2002. Regulation of noise in the expression of a single gene. Nat. Genet., 31:69–73.

28.

Paulsson

2004. Summing up the noise in gene networks. Nature, 427:415–418.

29.

Priestley

M.B.

1982. Spectral Analysis and Time Series. Academic Press: London.

30.

Proietti

F.A.

, Carneiro-Proietti

A.B.

, Catalan-Soares

B.C.

et al. 2005. Global epidemiology of HTLV-I infection and associated diseases. Oncogene, 24:6058–6068.

31.

Ramsey

, Orrell

, Bolouri

2005. Dizzy: stochastic simulation of large-scale genetic regulatory networks. J. Bioinform. Comput. Biol., 3:415–436.

32.

Rao

C.V.

, Wolf

D.M.

, Arkin

A.P.

2002. Control, exploitation and tolerance of intracellular noise. Nature, 420:231–237.

33.

Simpson

M.L.

, Cox

C.D.

, Sayler

G.S.

2003. Frequency domain analysis of noise in autoregulated gene circuits. Proc. Natl. Acad. Sci. USA, 100:4551–4556.

34.

Simpson

M.L.

, Cox

C.D.

, Sayler

G.S.

2004. Frequency domain chemical Langevin analysis of stochasticity in gene transcriptional regulation. J. Theor. Biol., 229:383–394.

35.

Smith

J.O.

2007. Mathematics of the Discrete Fourier TransformW3K Publishing [Online].http://books.w3k.org/. 2010 November 1.

36.

Taylor

G.P.

, Hall

S.E.

, Navarrete

et al. 1999. Effect of lamivudine on human T-cell leukemia virus type 1 (HTLV-1) DNA copy number, T-cell phenotype, and anti-Tax cytotoxic T-cell frequency in patients with HTLV-1-associated myelopathy. J. Virol., 73:10289–10295.

37.

Thattai

, van Oudenaarden

. 2001. Intrinsic noise in gene regulatory networks. Proc. Natl. Acad. Sci. USA, 98:8614–8619.

38.

Thomaseth

, Cobelli

1999. Generalized sensitivity functions in physiological system identification. Ann. Biomed. Eng., 27:607–616.

39.

T.L.

, Maheshri

2010. Noise can induce bimodality in positive transcriptional feedback loops without bistability. Science, 327:1142–1145.

40.

Weinberger

L.S.

, Dar

R.D.

, Simpson

M.L.

2008. Transient-mediated fate determination in a transcriptional circuit of HIV. Nat. Genet., 40:466–470.

41.

Wilkinson

D.J.

2006. Stochastic Modelling for Systems Biology. Chapman & Hall/CRC: London.

42.

, Xiao

, Ren

et al. 2006. Probing gene expression in live cells, one protein molecule at a time. Science, 311:1600–1603.

Sensitivity Analysis of Retrovirus HTLV-1 Transactivation

Abstract

Abstract

1. Introduction

2. Methods

2.1. Model of the HTLV-1-derived gene circuit

2.2. Modeling intrinsic noise

2.3. Sensitivity analysis

2.4. Gillespie's exact stochastic simulations

3. Results

4. Conclusion

5. Appendix

5.1. Molecular syntheses and degradations at steady state

5.2. HTLV-1 gene circuit transfer functions

5.3. Analytic calculation of ACFout

Footnotes

Acknowledgments

Disclosure Statement

References

5.3. Analytic calculation of ACF_out