Stochastic modeling of human papillomavirusearly promoter gene regulation

Abstract

High risk forms of human papillomaviruses (HPVs) promote cancerous lesions and are implicated in almost all cervical cancer. Of particular relevance to cancer progression is regulation of the early promoter that controls gene expression in the initial phases of infection and can eventually lead to pre-cancer progression. Our goal was to develop a stochastic model to investigate the control mechanisms that regulate gene expression from the HPV early promoter. Our model integrates modules that account for transcriptional, post-transcriptional, translational and post-translational regulation of E1 and E2 early genes to form a functioning gene regulatory network. Each module consists of a set of biochemical steps whose stochastic evolution is governed by a chemical Master Equation and can be simulated using the Gillespie algorithm. To investigate the role of noise in gene expression, we compared our stochastic simulations with solutions to ordinary differential equations for the mean behavior of the system that are valid under the conditions of large molecular abundances and quasi-equilibrium for fast reactions. The model produced results consistent with known HPV biology. Our simulation results suggest that stochasticity plays a pivotal role in determining the dynamics of HPV gene expression. In particular, the combination of positive and negative feedback regulation generates stochastic bursts of gene expression. Analysis of the model reveals that regulation at the promoter affects burst amplitude and frequency, whereas splicing is more specialized to regulate burst frequency. Our results also suggest that splicing enhancers are a significant source of stochasticity in pre-mRNA abundance and that the number of viruses infecting the host cell represents a third important source of stochasticity in gene expression.

Keywords: Feedback regulation; RNA splicing; Stochastic bursting.

Fig 1.. HPV molecular biology.

(A) HPV genome structure. The HPV genome is shown as a black circle with the early (P_E), late (P_L) and E8 (P_E8) promoters marked by arrows. The six early ORFs (in red), namely the regulatory genes E1, E2, E4, E5 are expressed from either P_E or P_L, the oncogenes E6 and E7 are expressed only by P_E and E8 is expressed only by P_E8, depending of the stage of the viral lifecycle. The late ORFs L1 and L2 (in orange) are expressed from P_L. All the viral genes are encoded on one strand of the double-stranded circular DNA genome. (B) The early promoter. This promoter accounts for positive and negative transcriptional regulation and is modeled with two binding sites where DE2 dimer and E1E2 heterodimer can bind and modulate transcription. (C) Splicing regulation. E1 and E2 transcript levels are spliced at the splicing sites SA742 and SA2709, respectively, while the remaining early genes are spliced at SA3358, regulated by the splicing factor SRSF1. E1 and E2 are produced in the absence of splicing factor SRSF1, hence they are indirectly regulated by this latter. (D) The SRSF1 promoter. This promoter is modeled as a two-state system with E2 transactivation. This constitutes a negative feedback loop for E1 and E2 expression. The formation of the heterodimer between SRSF1 protein and transcript, mSRSF1/SRSF1, is a negative regulator of SRSF1 activity.

Fig 2.. State diagram for the early promoter (P_E).

The promoter P_E,i is regulated by the binding of DE2 and E1E2. State P_E,0 (empty promoter) accounts for basal transcription, state P_E,1 (single bound DE2) accounts for positive feedback, state P_E,3 accounts for strong positive feedback and states P_E,2 (two bound DE2) and P_E,4 (E1E2 and DE2 bound) account for negative feedback. k_i, i ϵ {2,4,6,8} are dissociation rate constants while the k_i, i ϵ {1,3,5,7} are association rate constants.

Fig 3.. State diagram for the pM splicing.

pM is the pre-mRNA which is not undergoing splicing or it is undergoing splicing through splicing sites not accounted for in the model; pM1 is the pre-mRNA which is in a state where splicing occurs at SA742 site with the consequent conversion to E1 mRNA; pM2 is the pre-mRNA which is in a state where splicing occurs at SA2709 site with the consequent conversion to E2 mRNA; pM3 is the pre-mRNA which is in a state where splicing occurs at SA3358 site with the consequent conversion to the remaining early transcripts. The binding of SRSF1 splicing factor to the pM enhances the pM3 state occurrence. k_i, i ϵ {10,12,14} are dissociation rate constants while the k_i, i ϵ {9,11,13} are association rate constants.

Fig 4.. Steady-state probabilities for the early promoter.

Model predictions of the early promoter states as a function of the total copy number [CN] TE.

Fig 5.. Deterministic model behavior under nominal parameters.

Time series copy number (CN) for (A) pM primary transcript, (B) pM1 primary transcript, (C) pM2 primary transcript, (D) pM3 primary transcript, (E) mE1 transcript, (F) mE2 transcript, (G) E1 protein, (H) E2 protein, (I) DE2 homodimer, (J) E1E2 heterodimer.

Fig 6.. Features of the model response for pM1, pM2 and pM3 primary transcripts under different parameter values of the splicing control.

Features of the model response are investigated for different values of k₉ and k₁₃ parameters, which are the transition rates between pM and pM1 and pM and pM3, respectively. k₁₄, the transition rate between pM3 and pM, is also vaired being connected to k₁₃ through the constraint k₁₄ = k₁₃*k_D,SRSF1, where k_D,SRSF1 is the SRSF1 dissociation constant (see Table S1, Supplementary Information). k₁₃ axis, in all the subplots, is intended to be scaled by an order of magnitude of 10⁻⁷. (A) pM1 peak amplitude copy number [CN], (B) pM1 steady state [CN], (C) pM1 peak width [min], (D) pM1 peak time [min], (E) pM2 peak amplitude [CN], (F) pM2 steady state [CN], (G) pM2 peak width [min], (H) pM2 peak time [min], (I) pM3 peak amplitude [CN], (J) pM3 steady state [CN], (K) pM3 peak width [min], (L) pM3 peak time [min].

Fig 7.. Features of the model response for pM1, pM2 and pM3 primary transcripts under different parameter values of the splicing control and correspondent deterministic behavior.

Features of the model response are investigated for different values of k₁₀ and k₁₁ parameters, which are the transition rates between pM1 and pM and pM and pM2, respectively. (A) pM1 steady state copy number [CN], (B) pM2 steady state [CN], (C) pM3 steady state [CN]. Time series [CN] under a 10 fold decrease k₁₁ parameter for (D) pM primary transcript, (E) pM1 primary transcript, (F) pM2 primary transcript, (G) pM3 primary transcript, (H) mE1 transcript, (I) mE2 transcript, (J) E1 protein, (K) E2 protein, (L) DE2 homodimer, (M) E1E2 heterodimer.

Fig. 8.. Stochastic model behavior for slow promoter and splicing fluctuations regime.

Promoter and splicing binding rates were decreased of 10 fold compared to the nominal values. Time series of copy number, [CN], and steady-state distributions for: (A) pM primary transcript, (B) mE2 transcript, (C) DE2 homodimer, (D) E1E2 heterodimer. The comparisons between the stochastic behavior (red) and the quasi-equilibrium approximation (green) are shown during the transient response (left panels) and steady state (middle panels). Numerical steady state mean (blue) is also shown during the steady state. The right panels show the steady state distributions. Mean and distributions were derived from 100 stochastic realizations run for a time length of 4 days.

Fig. 9.. Stochastic model behavior for decreased SRSF1 promoter fluctuations.

SRSF1 promoter binding rates were decreased of 100 fold than the nominal value and the ratio between the SRSF1 synthesis rates between active and inactive promoter increased of 10 fold. Time series of copy number, [CN], and steady-state distributions for: (A) pM primary transcript, (B) mE2 transcript, (C) DE2 homodimer, (D) E1E2 heterodimer. The comparisons between the stochastic behavior (red) and the quasi-equilibrium approximation (green) are shown during the transient response (left panels) and steady state (middle panesl). Numerical steady state mean (blue) is also shown during the steady state. The right panels show the steady state distributions. Mean and distributions were derived from 100 stochastic realizations run for a time length of 4 days.

Fig. 10.. SRSF1 amplitude modulation.

Time series of copy number, [CN], and quasi-equilibrium for SRSF1 promoter fluctuations decreased of 100 fold than the nominal value and the ratio between the SRSF1 synthesis rates between active and inactive promoter increased of 10 fold. (A) pM primary transcript, (B) pM1 primary transcript, (C) pM2 primary transcript, (D) pM3 primary transcript, (E) DE2 homodimer, (F) SRSF1 protein. The stochastic behavior is shown in red and the quasi-equilibrium approximation is shown in blue.

Fig. 11.. Noise strength.

Coefficient of variation (CV) of some interesting chemical species under variation of the magnitude of the promoter binding rates (proportional to f_EP), splicing rates (proportional to f_SS) and SRSF1 synthesis rates (proportional to f_SRSF1). CV as a function of f_EP for: (A) pM1 primary transcript, (B) pM2 primary transcript, (C) mE1 transcript, (D) mE2 transcript. CV as a function of f_SS for: (E) pM primary transcript, (F) pM1 primary transcript, (G) pM2 primary transcript, (H) pM3 primary transcript. CV as a function of f_SRSF1 for: (I) pM primary transcript, (J) pM1 primary transcript, (K) pM2 primary transcript, (L) pM3 primary transcript.

Fig 12.. Promoter and splicing control of bursts amplitude and frequency.

Bursts features (frequency, duration and amplitude) of mE2 transcript are shown for different magnitudes of promoter binding rates (proportional to f_EP), splicing rates (proportional to f_SS), SRSF1 synthesis rates (proportional to f_SRSF1) and different ratios (R) between positive and negative feedback strengths. Bursts frequency (left panels), bursts duration (middle panels) and bursts amplitude (right panels). Bursts features as function of (A,B,C) f_EP, (D,E,F) R, (G,H,I) f_SS, (J,K,L) f_SRSF1,.

Fig 13.. Noise and bursts features controlled by the number of viruses N.

Coefficient of variation of (A) pM2 primary transcript and (B) mE2 transcript as function of N. (C-E) Burst frequency, duration and amplitude of mE2 transcript as a function of N.