The dynamics of the protein output can result in a faithful representation of the current biological environment.

<p>We consider a 2-stage model of gene expression <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi.1002965-Thattai1" target="_blank">[22]</a>. The extracellular environment or input, , gives the current rate of transcription and the signal of interest . We model as either a 2-state Markov chain with equal switching rates between states (the states each have unconditional probability of ) (A&C); or as proportional to a Poissonian birth-death process for a transcriptional activator (B&D; proportionality constant of 0.025). The transformed signals (in red, lower panels) are a perfect representation of , although protein levels (in blue) are not. , the lifetime of equals 1 hr, and the translation rate . Degradation rates of mRNA and protein are chosen to maximize the fidelity, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi.1002965.e163" target="_blank">Eq. 7</a>. The units for are chosen so that its variance equals one.</p

Dynamical error as the difference between two conditional expectations.

<p>To illustrate, we consider a 2-stage model of gene expression with the input, , equal to the current rate of transcription, and the signal of interest . We model as a 2-state Markov chain and show simulated trajectories of the protein output, , corresponding to four different input trajectories, . These input trajectories (or histories) all end at time in the state (not shown) and differ according to their times of entry into that state (labelled on the time axis; is off figure). (black lines) is the average value of at time given a particular history of the input : the random deviation of around this average is the mechanistic error (shown at time for the first realisation of ). is the average or mean value of given that the trajectory of ends in the state at time . (red line) can be obtained by averaging the values of over all histories of ending in . The mean is less than the mode of the distribution for because of the distribution's long tail. , not shown, is obtained analogously. The dynamical error, , is the difference between and and is shown here for the first trajectory, . <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi-1002965-g003" target="_blank">Fig. 3B</a> shows data from an identical simulation model (all rate parameters here as detailed in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi-1002965-g003" target="_blank">Fig. 3B</a>).</p

As the protein lifetime decreases, a trade-off between dynamical and mechanistic error determines fidelity.

<p>We consider a 2-stage model of gene expression with the input, , equal to the current rate of transcription, and the signal of interest . (A) The magnitude of the relative fidelity errors as a function of the protein degradation rate, (from <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi.1002965.e242" target="_blank">Eqs. 11</a>, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi.1002965.e247" target="_blank">12</a> and <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi.1002965.e248" target="_blank">13</a>), using a logarithmic axis. (B–D) Simulated data with as in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi-1002965-g001" target="_blank">Fig. 1A</a>. The units for are chosen so that its variance equals one in each case (hence and ). Pie charts show the fractions of the protein variance due to the mechanistic (m) and dynamical (d) errors and to the transformed signal. The latter equals . In B, the relative protein lifetime, , is higher than optimal () and fidelity is 2.2; in C, is optimal () and fidelity is 10.1; and in D, is lower than optimal () and fidelity is 5.3. Dynamical error, , is the difference between (black) and the faithfully transformed signal (red), and decreases from B to D, while mechanistic error increases. The lower row shows the magnitudes of the relative dynamical error (black) and relative mechanistic error (orange). All rate parameters are as in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002965#pcbi-1002965-g001" target="_blank">Fig. 1 A&C</a> with , unless otherwise stated.</p

Comparison of the Extrande and integral methods.

<p><b>(A)</b> Comparison of CPU times for Extrande and the modified next (MN) integral method [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.ref028" target="_blank">28</a>]. CPU times broken down into their constituents (color coded), and shown as a function of the look-ahead horizon, <i>L</i>, for Extrande (see also <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.box001" target="_blank">Box 1</a>). Time-step of input presimulation and of integration for the MN method both equal to 10⁻⁶h. CPU times were collected while simulating the two state model of gene expression with noisy circadian transcription (see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.g001" target="_blank">Fig 1D</a>) up to <i>t</i> = 10 days. <b>(B)</b> Percentage of exponential random variables generated in Step 4 of Extrande (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.box001" target="_blank">Box 1</a>) that are rejected, thinned, and accepted, as a function of <i>L</i>. Extrande simulation, network and input as in (A). <b>(C)</b> As in (A) but for the SynDM network (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.g003" target="_blank">Fig 3A</a>) with single OU input presimulated using a time-step of 10⁻²s (and with lifetime 1h, CV 0.5), and integration time-step for the MN method also equal to 10⁻²s. <b>(D)</b> Comparison of CPU times (for 10 simulated days) and of percentage errors for Extrande and the MN integral method. Network and input as in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.g001" target="_blank">Fig 1D</a> (and panels A & B), time-step of input presimulation again equal to 10⁻⁶h. The absolute value of the percentage error in the integral method’s estimate of the conditional mean is shown in red, both at 6h (crosses) and averaged over the first 24h (compared to Extrande, circles). CPU time for Extrande corresponds to an intermediate value of <i>L</i> (in practice, a few of the 1000 cells would be run initially to choose <i>L</i>). Throughout Fig 2, we use trapezoidal numerical integration for the MN integral method; the implementation of Extrande uses input presimulation over the look-ahead horizon <i>L</i> from which its ceiling value is obtained; and the CPU time for input presimulation is excluded since it is identical for the MN and Extrande methods.</p

The effect of extrinsic fluctuations from upstream quorum signaling on the competence decision of <i>B. subtilis</i>.

<p><b>(A)</b> The wild-type signaling (green) and competence (blue) modules. The Synthetic Decision-Making network (SynDM) has the additional positive regulation of ComK by pComA (dashed red arrow). Reaction networks and rate parameters described in detail in the <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.s002" target="_blank">S1 Text</a>. <b>(B–D)</b> Time courses of progress to competence shown for 100 cells containing the wild-type and SynDM networks, simulated using Extrande. In B & D, independent Gaussian, OU input processes for the pComA level in each cell are used, derived from an LNA model of the signaling module (see panel G). In C, pComA is held constant at the LNA mean of 1000 molecules. Progress to competence assumes differentiation proceeds with time-varying rate proportional to the level of ComK (see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.s002" target="_blank">S1 Text</a>), with progress equal to 1 corresponding to entry to competence. At time zero, the level of ComS and ComK mRNAs and proteins set to zero. <b>(E)</b> For the wild-type and SynDM networks, the percentage change in the fraction of a population of 1000 quorum sensing cells entering competence (within 20 hours) compared to the fraction when pComA is constant at 1000 molecules, as a function of the lifetime and CV of the OU input modeling pComA fluctuations. The limit with pComA constant in each cell is also shown, drawn from a Gaussian distribution with mean 1000 and the indicated CV. <b>(F)</b> For the wild-type and SynDM networks, the estimated Prob[Competence|〈pComA〉] as a function of 〈pComA〉, the time-averaged level of pComA over the 20h experiment, for different OU inputs modeling pComA fluctuations. Estimation performed using logistic regression. <b>(G)</b> For the wild-type and SynDM networks, the fraction of a population of 1000 quorum sensing cells entering competence as a function of the proportionality constant of ComK-driven differentiation, for different OU inputs modeling pComA fluctuations (lifetime of 28s corresponds to model of upstream signaling lacking gene expression of the component proteins). <b>(H)</b> The autocorrelation function of pComA given by the LNA model of upstream signaling, compared to that of a single OU input process and 2 summed, independent OU processes, both having the mean and variance of pComA given by the LNA.</p

Comparison of the accuracy of the Extrande and SIA simulation methods.

<p><b>(A)</b><b>Gene expression with circadian transcription rate</b>, <i>k</i>(<i>t</i>)/<i>k</i>_<i>dp</i> = 4(1 + sin(2<i>πft</i>)), and period <i>f</i>⁻¹ = 24h. The proportional root mean square error between the average protein number from the SIA method, 〈<i>n</i>(<i>t</i>)〉, and the exact, time-dependent solution, <i>n</i>_ex(<i>t</i>), is shown as a function of the relative frequency of oscillation. The error is given by , where is the time-average of the exact solution. Physiological parameters for circadian rhythms (CR) in 4 different cell types are indicated. Actual mean protein numbers are varied via the translation rate (<i>k</i>_<i>s</i>), holding degradation rates constant. The error is particularly conspicuous (>60%) for <i>O. tauri</i> over the whole range of average protein numbers (10–10,000) whereas the exact Extrande method (Inset) accurately predicts the mean protein numbers for this case within sampling error (given by the standard error of the mean, SEM). <b>(B)</b> <b>Gene expression with noisy transcription rate</b>, <i>k</i>(<i>t</i>) = 〈exp(<i>ξ</i>(<i>t</i>)〉⁻¹exp(<i>ξ</i>(<i>t</i>)) where <i>ξ</i>(<i>t</i>) is zero-mean Gaussian (OU) noise with autocovariance 〈<i>ξ</i>(<i>t</i>)<i>ξ</i>(<i>t</i>′)〉 = 5<i>e</i>^{−<i>γ</i>|<i>t</i>−<i>t</i>′|}. We show the proportional error of the stationary average protein number from the SIA method as <i>γ</i> varies, with average protein numbers set via <i>k</i>_<i>s</i> as in (A). Autocorrelation times of the transcription rate of the order of the cell cycle (CC) are indicated for 4 different cell types. The error is particularly conspicuous (60–90%) for stable proteins removed mainly by dilution, as is common in bacteria (<i>γ</i>/<i>k</i>_dp = 1), where we show (Inset) simulated average protein numbers for a population of 100 cells. The error bars denote one standard deviation of the bootstrap distribution. <b>(C)</b> <b>Noisy circadian oscillations in an <i>O</i>. <i>tauri</i> cell population</b>. Average protein numbers for 2500 cells and 10 days simulated using a circadian transcription rate with cell cycle-induced amplitude fluctuations on a similar timescale: <i>k</i>(<i>t</i>) = 20exp(<i>ξ</i>(<i>t</i>))(1 + sin(2<i>πft</i>)), where <i>ξ</i>(<i>t</i>) is zero-mean Gaussian (OU) noise with autocovariance and <i>γ</i>/<i>k</i>_dp = <i>f</i> ln 2. While Extrande correctly predicts sustained oscillations (blue), the SIA method predicts only damped oscillations (red). Extrande is in excellent agreement with the corresponding moment equations of the master equation (dots, equivalent to ODE solution). Single cell realizations (Inset) reveal the SIA method shows unphysical loss and revival of oscillations. <b>(D)</b> <b>Average behaviour of <i>O</i>. <i>tauri</i> cells conditional on transcription dynamics</b>: We pregenerated a single realization (Inset) of the transcription rate, <i>k</i>(<i>t</i>), used in (C), and averaged over 1,000 resultant protein trajectories (all parameters as in C). The solution of the corresponding SDE for average protein conditional on the trajectory of <i>k</i>(<i>t</i>) agrees very well with the average from Extrande, in contrast to the SIA method. See <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004923#pcbi.1004923.s002" target="_blank">S1 Text</a> for simulation details and other rate parameters.</p

Supplementary figures.

Fig A: Testing HormoneBayes on synthetic data. Fig B: Assessing the effect of the prior for the LH clearance rate. Fig C: Tuning HormoneBayes when pulses are not clear by using a more informative prior on parameter f. Fig D: Pulse identification using HormoneBayes. Fig E: Using HormoneBayes to identify the effect of interventions on LH pulsatility. (PDF)</p

pcbi.1011928.t001 - HormoneBayes: A novel Bayesian framework for the analysis of pulsatile hormone dynamics

pcbi.1011928.t001 - HormoneBayes: A novel Bayesian framework for the analysis of pulsatile hormone dynamics</p

HormoneBayes handles LH pulsatility analysis in different contexts.

(A) Inferred pulsatility strength and maximum secretion rate parameters for different individuals: healthy men (n = 10); healthy post-menopause women (n = 13); healthy pre-menopausal women (n = 4). (B) Inferred parameters for healthy pre-menopausal women (n = 4); women with PCOS (n = 6) and women with HA (n = 5) illustrating how the assessment of LH pulsatility could help facilitate diagnosis of patients presenting with reproductive endocrine disorders. (C) Representative fits of the model are given for one subject in each dataset.</p

Pulse identification using HormoneBayes.

(A) Pulses can be identified using the expected value of the pulsatile hypothalamic signal, which can be interpreted as the probability of a pulse at a given timepoint. Here, we mark the onset of a pulse when the pulsatile hypothalamic signal crosses the 0.5 threshold, indicating that at this point a pulse is the most likely event. (B) The majority of the identified pulses (89%, 77/87) are in line with those obtained using the deconvolution method. For the analysis we used LH data obtained from healthy pre-menopausal women in early follicular phase (n = 16).</p