Global sensitivity analysis of a dynamic model for gene expression in Drosophila embryos

Gene regulation in embryonic development

Embryonic development in animals is very precisely controlled by a network of regulatory proteins (Davidson, 2010; Peter and Davidson, 2011). For any particular protein, the exact level of expression at a specific time point can be crucial to the proper development of the organism (Davidson & Levine, 2008). Within each cell of the developing embryo, protein abundance is a function of two key molecular events: transcription and translation. Transcription is the process of reading a gene in a DNA template to produce a messenger RNA (mRNA), while translation is the process of reading the mRNA to produce a protein product. A simple, but long-standing question in Biology is the following: Which is contributing more to the variance in protein levels in cells, transcription or translation?

Work analyzing the importance of transcription in mammalian cells

Many experimental studies have been conducted in an attempt to understand the roles of transcription and translation in regulating the dynamic nature of mRNA and protein concentrations (Maier, Guell & Serrano, 2009; Vogel et al., 2010; Schwanhausser et al., 2011; Beck et al., 2011; Bantscheff et al., 2012; Vogel & Marcotte, 2012; Li, Bickel & Biggin, 2014). One such recent detailed study, conducted by the Biggin lab using data from ubiquitously expressed (housekeeping) genes in cultured mammalian cells, aimed to improve existing quantification of protein abundances through statistical analysis of the impact of experimental error (Li, Bickel & Biggin, 2014). In this study, a two-part regression procedure was used to derive new estimates of protein abundance from the 2011 data set of Schwanhausser et al. (2011). Using these new, corrected measurements of protein abundances, along with previous measurements of mRNA and protein degradation rates, they were able to determine the relative importance of transcription, translation, mRNA degradation, and protein degradation (see Fig. 1A). This analysis is referred to as the measured protein error strategy. The result of this procedure found that for the 4,212 genes considered, transcription contributed the most to protein abundance (∼38%), with translation contributing slightly less (∼30%), followed by mRNA degradation (∼18%) and protein degradation (∼14%) (Li, Bickel & Biggin, 2014). This result is important to note, as it differs drastically from Schwanhausser et al.’s original conclusion that translation accounts for the largest contribution (∼55%) to overall variance in the cellular abundance of proteins (Schwanhausser et al., 2011; Li, Bickel & Biggin, 2014).

Existing mathematical models of gene expression

Due to the quantitative nature of gene regulation in the embryo, and the advent of new experimental techniques giving rise to massive amounts of mRNA and protein concentration data, various mathematical models have been derived and implemented to help understand the complexity that lies within developmental gene regulatory networks. These models range from static models, considering only transcription at a single time point in development in a single cell, to dynamic spatio-temporal models that incorporate transcription, translation, diffusion, and decay rates for a network of genes that regulate one another over a continuous time frame (Jaeger et al., 2004; Santillan & Mackey, 2004; Bintu et al., 2005; Janssens et al., 2006; Zinzen et al., 2006; Segal et al., 2008; Gertz, Siggia & Cohen, 2009; Fakhouri et al., 2010; He et al., 2010; Bieler, Pozzorini & Naef, 2011; Janssens et al., 2013; Ilsley et al., 2013; Dresch et al., 2013; Samee & Sinha, 2014).

To accurately model protein abundance in a metazoan animal, such as a mouse or fruit fly, one must consider both spatial and temporal dynamics in the developmental system. One such model that we developed uses a discretized reaction–diffusion equation to model concentrations of mRNA and protein in a developing Drosophila embryo across n nuclei (Dresch et al., 2013). This model not only incorporates terms for mRNA and protein synthesis and decay, but also diffusion of these molecules in the developing embryo. This is particularly important in Drosophila development as the early stages of embryogenesis are marked by 13 mitotic (nuclear) divisions in the absence of cellular divisions, resulting in a multinucleate syncytial embryo (see Fig. 1B).

In its simplest form, the model can be written as: (1)${\displaystyle \frac{d{y}_{a,i}}{dt}=S_{a,i}\left(Y\right)+D_a\left(y_{a,i+1}-2y_{a,i}+y_{a,i-1}\right)-{\lambda }_a{y}_{a,i}}$ where a represents the specific mRNA or protein that the equation corresponds to, i represents the specific nucleus, D_a and λ_a are the corresponding diffusion and decay rates, and Y represents the entire vector of mRNA and protein concentrations within the system being modeled.

Many similar models have been used to model the expression dynamics of the gap gene system in the developing Drosophila embryo (Jaeger et al., 2004; Okabe-Oho et al., 2009; Ashyraliyev et al., 2009; Bieler, Pozzorini & Naef, 2011; Holloway et al., 2011; Janssens et al., 2013; Holloway & Spirov, 2015). Although these models all rely on an underlying reaction–diffusion framework, they vary greatly in their implementation. Both deterministic (Jaeger et al., 2004; Ashyraliyev et al., 2009; Bieler, Pozzorini & Naef, 2011; Janssens et al., 2013) and stochastic (Okabe-Oho et al., 2009; Holloway et al., 2011; Holloway & Spirov, 2015) models have been able to accurately predict the effects of particular perturbations to the network. Stochastic models of hunchback regulation have been used to shed light on the underlying factors that reduce noise and promote stability of the hunchback gradient. These include the number or arrangement of BICOID and KRUPPEL binding sites within the regulatory sequences that control transcription of the hunchback gene as well as protein diffusion (Okabe-Oho et al., 2009; Holloway et al., 2011; Holloway & Spirov, 2015). In this study, we focus on the broad impacts of transcription, translation, diffusion, and decay, and do not consider specific transcription factor binding sites within regulatory DNA sequences. Thus, we focus on a deterministic model and do not include any stochasticity.

Global parameter sensitivity analysis and HDMR

To help develop a better understanding of the model parameters, including how their values impact the model output (protein abundance) and how one might interpret that impact with respect to the biological system, parameter sensitivity analysis is needed. Parameter sensitivity analysis refers to a mathematical analysis of the change in model output as a result of variation in the input parameter values (Frey & Patil, 2002; van Riel, 2006; Tang et al., 2006; Dresch et al., 2010; Ay & Arnosti, 2011; Jarrett et al., 2014). This analysis can be done locally, at a particular point in parameter space, or globally across the entire parameter space.

Local parameter sensitivity analyses are typically implemented by simply computing or approximating the partial derivative of the objective function at a particular point in parameter space to determine how the function changes locally with respect to small variations of a particular parameter (van Riel, 2006; Dresch et al., 2010; Reeves & Fraser, 2009). The major advantages to adopting a local approach are that it is straightforward, often easy to interpret, and computationally inexpensive. However, a significant limitation of local methods is that when dealing with a large parameter space focusing on a single point in that space may not be representative of much of the overall parameter space. In contrast, global methods allow one to calculate parameter sensitivities while considering the full range of parameter space (Jarrett et al., 2014). Most global methods also have the ability to calculate higher order sensitivities, which capture interactions between parameters. This can be challenging to do with a local method, especially when parameter space is large and many different parameter combinations within that space lead to valid model outputs. Therefore, in this study we focus on a global method for parameter sensitivity analysis for our model of gene expression in the Drosophila embryo.

Higher Dimensional Model Representation (HDMR) is a robust global method for calculating parameter sensitivities (Ziehn & Tomlin, 2009; Dresch et al., 2010). This method uses the Monte Carlo integration method to decompose the model output, f(x), often referred to as the objective function, into terms of increasing dimensionality with respect to the parameters x₁, …, x_N: (2)${\displaystyle f\left(x\right)=f_0+\sum _{i=1}^N{f}_i\left(x_i\right)+\sum _{1\le i<j\le N}{f}_{ij}\left(x_i,x_j\right)+\cdots +f_{12\dots N}\left(x_1,x_2,\dots ,x_N\right).}$

In the above equation, f₀ is the main effect, and is approximated by the overall mean of model output over all parameter sets sampled. Each function f_i(x_i) is a first-order term representing the effect of the parameter x_i acting independently on the model output. Each function f_ij(x_i, x_j) is a second-order term representing the effect of the parameters x_i and x_j on the model output. These terms represent the impact pairwise parameter interactions have on determining the model output.

This approximation is done on a bounded subset of ℝ^N, where N is the number of model parameters. The bounded subset represents the parameter space and in each dimension, corresponds to a realistic range for the given parameter. In some applications, the parameter space is known or experimentally determined using empirically determined measurements (Ziehn & Tomlin, 2009; Tang et al., 2006). In other cases, parameter space is chosen based on model assumptions and simulations (Frey & Patil, 2002; Gutenkunst et al., 2007; Dresch et al., 2010; Ay & Arnosti, 2011; Jarrett et al., 2014).

One of the main assumptions when using HDMR is that the objective function values are normally distributed (Ziehn & Tomlin, 2009). Thus, using the bounded parameter space, one generates a pseudo-random sampling using a Sobol Set (Sobol, 1976). The bounded parameter space must be a hypercube in N-dimensional space, so it can be normalized to the unit hypercube for sampling. Each set of parameter values is then used as input for a model simulation and the corresponding model output is obtained. Once this sampling has been done for all parameter sets sampled, HDMR approximates the mean of the output values, f₀, as well as higher order terms using orthonormal polynomial approximations. First- and second-order terms are then normalized by the total variance to obtain the main effect of each parameter and the effects of pair-wise parameter interactions, referred to as first- and second-order sensitivity indices respectively. Although this method has the ability to calculate higher-order sensitivities, it has been shown that using only first- and second-order terms is sufficient to approximate the total sensitivity in the system (Li, Rosenthal & Rabitz, 2001; Liang & Guo, 2003; Dresch et al., 2010).

In this study, we utilize a global HDMR analysis to investigate the sensitivity of our Drosophila embryo gene expression model to the individual transcription, translation, diffusion, and decay rate parameters and higher-order interactions between these parameters. In addition, we compare our results to those in mammalian cells and other studies that have attempted to model different gene expression systems (Li, Bickel & Biggin, 2014).

Simplified model and parameters

In this study, we use a simplified version of our earlier model (Dresch et al., 2013), which was used to predict both mRNA and protein concentrations along a one-dimensional strip of nuclei in a developing Drosophila embryo.

The simplifying assumption applied in the current study is that the gene of interest has spatially uniform transcriptional activity. Thus, the transcription rate is held constant. This allows us to utilize simple numerical solvers such as Euler’s method or Runga-Kutta methods, and to measure the relative importance of transcription to the model output using a single parameter, σ. Note that diffusion is discretized with respect to space and zero flux boundary conditions are used (Dresch et al., 2013). For 2 ≤ i ≤ n − 1, where n corresponds to the number of nuclei being modeled, the model equations are thus as follows: ${\displaystyle \frac{d{y}_1}{dt}=\sigma +d\left(y_2-y_1\right)-\lambda y_1}$ ${\displaystyle \frac{d{y}_i}{dt}=\sigma +d\left(\left(y_{i+1}-y_i\right)+\left(y_{i-1}-y_i\right)\right)-\lambda y_i}$ ${\displaystyle \frac{d{y}_n}{dt}=\sigma +d\left(y_{n-1}-y_n\right)-\lambda y_n}$ ${\displaystyle \frac{d{y}_{n+1}}{dt}=\tau y_1+\delta \left(y_{n+2}-y_{n+1}\right)-\gamma y_{n+1}}$ ${\displaystyle \frac{d{y}_{n+i}}{dt}=\tau y_i+\delta \left(\left(y_{n+i+1}-y_{n+i}\right)+\left(y_{n+i-1}-y_{n+i}\right)\right)-\gamma y_{n+i}}$ ${\displaystyle \frac{d{y}_{2n}}{dt}=\tau y_n+\delta \left(y_{2n-1}-y_{2n}\right)-\gamma y_{2n}.}$

Here, y_j represents mRNA concentrations for 1 ≤ j ≤ n and protein concentrations for n + 1 ≤ j ≤ 2n. A schematic of the biological processes incorporated in the model is shown in Fig. 1. The model parameters are defined in Table 1. For all of the analysis and results shown in Table 2 and Figs. 2–4, 52 nuclear positions are used at approximately 50% egg height (ventral-dorsal) from approximately 10% to 90% egg length (anterior-posterior) and the numerical solver Runga-Kutta 4 is used to approximate the solutions to the system of Ordinary Differential Equations shown above.

Exploring parameter space

Before a global parameter sensitivity analysis can be performed, one must first define parameter space. This is done by finding a range for each model parameter that results in ‘realistic’ model predictions. This will result in a six-dimensional hypercube, as required for the Sobol Set sampling method (Sobol, 1976). For our model, we determine this range for each parameter by exploring six-dimensional parameter space and recording the parameter value combinations resulting in model predictions of reasonable protein concentrations.

The exploration of parameter space is done iteratively on a parameter-by-parameter basis. First, we define all protein concentration values ≥7 as protein saturation and ≤0.01 as undetectable protein. The method then works in the following way:

1. Initial parameter ranges are set to $\sigma ,d,\lambda ,\delta \in \left[0,1\right]$, $\tau \in \left[0.012,1\right]$, and $\gamma \in \left[0.2248,1\right]$. Note: Lower bounds were chosen using the values estimated in Dresch et al. (2013).

2. Choose a single parameter. Test each of the boundaries of the parameter range by holding the parameter of interest constant at the lower (or upper) bound and varying all other parameters within their current ranges. Then, run all simulations for $t\in \left[0,10\right]$ with all 5 different initial conditions. If >5% of all simulations result in saturated or undetectable protein concentrations, then modify the boundary of the parameter range by increasing or decreasing it by a number between 0.001 and 0.1.

   Note: Lower bounds on parameter ranges were never allowed to go below 0.0, as negative parameter values would violate the model assumptions.

3. Repeat step 2 for all remaining parameters.

4. When all parameter ranges have been tested, let all 6 parameters vary within their current ranges. If >5% of all simulations result in saturated or undetectable protein concentrations, then go back to step 2. If not, then stop.

The realistic parameter ranges we obtained and used for our subsequent sensitivity analyses are listed in Table 1.

Objective function

To calculate model parameter sensitivies, we utilize a previously developed HDMR MATLAB script (Ziehn & Tomlin, 2009). Due to the spatial and temporal dynamics of our model, we perform this analysis using a variety of different objective functions.

We focus our analysis on a twenty minute time window in an early blastoderm embryo. Within this time window, we perform our sensitivity analysis at eleven distinct, evenly distributed time points. This removes any bias in the time point chosen and allows us to look at the behavior of parameter sensitivities over time.

One should note that a parameter’s first-order sensitivity is calculated by approximating the variance with respect to that parameter divided by the total variance of the objective function. Thus, at t = 0, this results in a ratio of $\frac{0}{0}$ since all parameter sets will result in an output equal to the initial condition. Thus, we define the first-order sensitivity of each parameter to be zero when the model output used is from the time point t = 0.

Spatially, the model aproximates mRNA and protein concentrations at 52 nuclei across the anterio-posterior axis of the embryo. When performing our sensitivity analysis, we focus on protein concentrations only and consider three different spatial locations: a nucleus in the anterior of the embryo, in the posterior of the embryo, or in the center of the embryo. We also perform sensitivity analyses using the mean protein concentration over all 52 nuclei.

To avoid any inherent bias in our results, we compute parameter sensitivities using each of the five initial conditions and each combination of spatial and temporal protein concentrations. Thus, the results shown are represenatitive of 120 runs of the HDMR algorithm.