Communicating oscillatory networks: frequency domain analysis

Background Constructing predictive dynamic models of interacting signalling networks remains one of the great challenges facing systems biology. While detailed dynamical data exists about individual pathways, the task of combining such data without further lengthy experimentation is highly nontrivial. The communicating links between pathways, implicitly assumed to be unimportant and thus excluded, are precisely what become important in the larger system and must be reinstated. To maintain the delicate phase relationships between signals, signalling networks demand accurate dynamical parameters, but parameters optimised in isolation and under varying conditions are unlikely to remain optimal when combined. The computational burden of estimating parameters increases exponentially with increasing system size, so it is crucial to find precise and efficient ways of measuring the behaviour of systems, in order to re-use existing work. Results Motivated by the above, we present a new frequency domain-based systematic analysis technique that attempts to address the challenge of network assembly by defining a rigorous means to quantify the behaviour of stochastic systems. As our focus we construct a novel coupled oscillatory model of p53, NF-kB and the mammalian cell cycle, based on recent experimentally verified mathematical models. Informed by online databases of protein networks and interactions, we distilled their key elements into simplified models containing the most significant parts. Having coupled these systems, we constructed stochastic models for use in our frequency domain analysis. We used our new technique to investigate the crosstalk between the components of our model and measure the efficacy of certain network-based heuristic measures. Conclusions We find that the interactions between the networks we study are highly complex and not intuitive: (i) points of maximum perturbation do not necessarily correspond to points of maximum proximity to influence; (ii) increased coupling strength does not necessarily increase perturbation; (iii) different perturbations do not necessarily sum and (iv) overall, susceptibility to perturbation is amplitude and frequency dependent and cannot easily be predicted by heuristic measures. Our methodology is particularly relevant for oscillatory systems, though not limited to these, and is most revealing when applied to the results of stochastic simulation. The technique is able to characterise precisely the distance in behaviour between different models, different systems and different parts within the same system. It can also measure the difference between different simulation algorithms used on the same system and can be used to inform the choice of dynamic parameters. By measuring crosstalk between subsystems it can also indicate mechanisms by which such systems may be controlled in experiments and therapeutics. We have thus found our technique of frequency domain analysis to be a valuable benchmark systems-biological tool.


Supplementary results
Crosstalk between oscillatory models of the cell cycle G1/S phase, NF-κB and p53

Increased coupling strength
The coupling strength of the oscillatory systems is controlled by the values of rate constants R30b and R40b. The initial coupling strength was chosen to be plausible and to respect known parameters and experimentally verified behaviour. To validate our results, we also considered twice and ten times this coupling strength. Figure S1 illustrates the perturbation (D) of cell cycle components by p53a using double coupling strength. While the trend of crosstalk is reinforced by most components, in line with intuition, it is noticeable that E2F-Rbpp is apparently influenced less with increased coupling strength (D = 0.178 vs. 0.476). With ten times coupling strength ( Figure S4) this unexpected trend is partially reversed (0.369). This counter-intuitive behaviour is to be expected in non-linear dynamical systems that are sensitive to both the amplitude and frequency of perturbations. Overall, Figures S1-S6, variously illustrating the influence of NF-B alone, p53 alone and the fully coupled systems with double and ten times coupling strength, tend to follow a more or less intuitive trend of increased measured perturbation with increased coupling. Whereas with single coupling strength the influence of p53a tends to be local and the influence of NF-Bn tends to be indirect, with increased coupling strength more species are affected to a greater degree and this demarcation is blurred. There are some noticeable exceptions, however. Species E2F-Rb, p16, Rb and Skp2 are minimally perturbed by any combination of coupling and coupling strengths.
Values for double and ten times coupling strength are given in Tables S2 and S3, respectively.  Table S2. Values are given in Table S2. Values are given in Table S2.  Figure S4: Perturbation of the cell cycle by p53a with ten times coupling strength.
Values are given in Table S3. Figure S5: Perturbation of the cell cycle by NF-Bn with ten times coupling strength.
Values are given in Table S3. Values are given in Table S3.   Table S3: Crosstalk in the cell cycle coupled to p53, to NF-κB and both p53 and NF-κB with ten times coupling strength. Each species of the cell cycle is assigned a D value (Equation (4)

Supplementary methods
The following subsections describe the construction of the system of networks and simulation models used in the main text. The three networks are given as independent models, with the coupling reactions itemised separately. The rate constants and initial numbers of molecules are common to both the stochastic and quasi-deterministic models. The following general assumptions were made: earlier studies revealed oscillatory NF-B activity in cells lacking IB  &  [1] and biphasic dynamics in cells in which NF-B-inducible IB was over-expressed [2]; induced expression of IB also gives rise to an oscillatory NF-B signal that is out of phase with IBinduced oscillations [3] and which helps to keep the late phase of TNF-and IL1--induced NF-B activity steady [4]; possibility that oscillations are a trade-off against rapid response to inflammatory signals and the necessity of additional feedback to provide steady supply of active NF-B, which might arguably be present in other linking pathwayse.g. p53-MDM2; in contrast to some other models that focus mainly on G1/S-phase transition control [5][6][7][8][9], this model is only interested in key processes involving data sets previously published (which strengthens the idea that cell cycle process can also be induced by NF-B active proteins).

NF-κB system
We and others have previously re-constructed a computational model to explain NF-B activation events following IKK activation by TNF stimulation [2,10]. catalytic degradation of the NF-B Inhibitor proteins due to IKK-induced phosphorylation and subsequent ubiquitination.

Cell cycle system
Many groups have reported the construction of the mammalian cell cycle models, the most recent being [11,12].

DNA damage transduction system
Stimuli such as DNA damage can activate both the NF-κB and the p53 pathways [13]. While p53 induces cell-cycle arrest or cell death in response to these treatments, the contribution of NF-κB to cell fate is more complex, and pathways in which it either antagonizes or cooperates with p53 have been described. NF-κB-mediated negative regulation of p53 can contribute to tumorigenesis and has been shown to operate at a number of levels. The delay oscillator describe by Geva-Zatorsky et. al [14] was the model of choice for the p53 system.

Stochastic delay differential equations
Unlike the cell cycle and NF-κB models, which are described by standard chemical and enzymatic reactions that can be straight-forwardly converted into the stochastic domain, the original p53 model is constructed around delay differential equations (DDEs) and requires special consideration when simulating stochastically. DDEs allow the current rate of a reaction to be dependent on the state of the system at some time in the past, abstracting potentially very complex (unknown) mechanisms into a combination of delays. Stochastic simulations based on a variant of the Gillespie algorithm [15] treat the evolution of the system as a Markov process, such that the current rate of a reaction is only dependent on the current state of the system. To incorporate a model based on DDEs into a stochastic simulation of this kind, it is also necessary to store the past states of the system. We thus implemented a function in our simulation software which returns the number of molecules of species X at τ time units in the past: delay(X,τ). Given that the time course of a stochastic simulation consists of the numbers of molecules of different species recorded at irregular time points 0, t 1 , t 2 , ... etc., where 0 < t 1 < t 2 < ..., the amount of species X at time t can be described by the sequence X 0 , Xt 1 , Xt 2 , ..., etc. At time t, the value returned by delay(X,τ) is then Xt i , where i is the maximum value that satisfies t i ≤ t-τ. If t < τ (as may happen at the beginning of a simulation), delay(X,τ) returns X 0 . This algorithm is consistent with the standard deterministic interpretation of delay differential equations and guarantees that for any specified initial state and past state corresponding to the delay used, the magnitude and direction of the average rate of leaving the state in the stochastic and quasi-deterministic models is identical (allowing for the conversion from concentration to numbers of molecules) to that for the deterministic case.
For the purposes of simulation we used simplified names of the chemical species. The following table maps the names used in the text to the names used in the models.

Stochastic (reaction-based) models
The following models are based on standard chemical reactions of the kind A + B  C + D, where a single reaction event simultaneously consumes molecules of A and B while producing molecules of C and D. Where not otherwise stated, the rates of reactions are calculated using the given constants and the assumption of mass action. An explicit rate function to generate the reaction propensity [15] is given for reactions where this does not apply. The function delay(·) is defined above. The symbol Ø is used to denote an arbitrary source or sink in creation and degradation reactions, respectively.

Quasi-deterministic (combined production / consumption) models
The following models are based on quasi-deterministic 'reactions', where a single reaction event is an unsynchronised production or consumption of a single molecule of a given species. The propensity [15] of a quasi-deterministic reaction is given by an explicit rate function that is the signed sum of the propensities of creation and consumption reactions.

Rate and other constants
The following constants are common to both sets of models. The reaction rates are derived from models based on concentration, hence alpha has units of l mol -1 and is the constant which is used to

Initial numbers of molecules
The following table contains the initial numbers of molecules used by all the simulation models.
Note that the values are either 0 or some number (an initial concentration in units of mol l -1 ) multiplied by alpha, the global constant used to specify the number of molecules in the system.

Supplementary example: a stochastic model of the eukaryotic cell cycle
The presented technique of frequency domain analysis can be particularly useful and revealing when applied to stochastic simulations of chemical systems containing one or more species in low copy numbers. The apparent behaviour in such simulation time courses may be very noisy, yet the underlying average behaviour will nevertheless be obvious to human observers. Frequency domain analysis is a means to formalise the perceived behaviour. Figure S7 illustrates the differences between stochastic and deterministic simulations, using as example the generic model of the eukaryotic cell cycle in [18]. Table S5 describes the stochastic model, extracted from the ODEs, comprising elemental reactions using mass action kinetics, enzymatic reactions with arbitrary kinetics and mass balance equations. Parameters are chosen to represent budding yeast and the species names are those used in [18]. The deterministic model exhibits limit cycle oscillation, making it conceivable to run arbitrarily long stochastic simulations and so arbitrarily define the resolution of the frequency domain analysis. This is not in general guaranteed: when a continuous deterministic model is discretised and made stochastic (or quasi-deterministic, as described above), it may contain states which have a non-zero probability of being reached but from which the system cannot exit. These absorbing states may exist in reality or may be unforeseen artefacts of the ODE approximation of reality, hence the validity of stochastic models created from deterministic systems containing arbitrary simplifications and abstractions is sometimes questioned. Such questions may be answered by the presented methodology.  fully stocahastic (red) models of the budding yeast cell cycle. The deterministic spectrum (created from a single time course of 10000 minutes sampled at 5 minute intervals) is clearly 'spiky' in nature, with many evident high frequency components and apparent numerical artefacts. By contrast, the average stochastic frequency spectrum (red) contains only four discernable low frequency peaks that are relatively rounded. The spectrum of the quasi-deterministic simulation appears closer to the fully stochastic than to the deterministic, however it contains three more discernable peaks and at higher frequencies it follows more closely the trend of the deterministic spectrum. The peaks of the stochastic and quasi-deterministic spectra apparently align with peaks in the deterministic spectrum, suggesting that the three systems have the same average primary mode of oscillation, however this alignment between models is not in general guaranteed. In each case the spectral value at zero frequency corresponds to the long term mean of the time series.  Table S5. As a result of random phase shifts between simulation runs, oscillatory behaviour in the average trace decays with time and hovers around the long term average number of molecules of CycBT (grey line). B Average frequency distribution of CycBT in stochastic (red), quasi-deterministic (blue) and determinist (black) models. The deterministic spectrum is clearly more 'spiky' than the stochastic spectra and has many high frequency components that are apparently 'lost in the noise' of the other models. The average time course was created from 800 simulation traces of 4000 minutes sampled at 2 minute intervals. The average frequency spectra were created from 800 traces of 10000 minutes sampled at 5 minute intervals. The deterministic spectrum was created from a single trace of 10000 minutes sampled at 5 minute intervals.    [18] was resolved into elemental reactions with mass action kinetics and enzymatic reactions having arbitrary kinetic laws. A constant of alpha = 424 l mol -1 was used to convert initial concentrations and rate constants to numbers of molecules. To discretise the cell growth, a mass granularity constant of beta = 1000 was adopted.