A model invalidationbased approach for elucidating biological signalling pathways, applied to the chemotaxis pathway in R. sphaeroides
 Mark AJ Roberts^{1, 2, 3},
 Elias August^{1, 3},
 Abdullah Hamadeh^{1, 3},
 Philip K Maini^{3, 4},
 Patrick E McSharry^{3, 5},
 Judith P Armitage^{2, 3} and
 Antonis Papachristodoulou^{1, 3}Email author
DOI: 10.1186/175205093105
© Roberts et al; licensee BioMed Central Ltd. 2009
Received: 02 June 2009
Accepted: 31 October 2009
Published: 31 October 2009
Abstract
Background
Developing methods for understanding the connectivity of signalling pathways is a major challenge in biological research. For this purpose, mathematical models are routinely developed based on experimental observations, which also allow the prediction of the system behaviour under different experimental conditions. Often, however, the same experimental data can be represented by several competing network models.
Results
In this paper, we developed a novel mathematical model/experiment design cycle to help determine the probable network connectivity by iteratively invalidating models corresponding to competing signalling pathways. To do this, we systematically design experiments in silico that discriminate best between models of the competing signalling pathways. The method determines the inputs and parameter perturbations that will differentiate best between model outputs, corresponding to what can be measured/observed experimentally. We applied our method to the unknown connectivities in the chemotaxis pathway of the bacterium Rhodobacter sphaeroides. We first developed several models of R. sphaeroides chemotaxis corresponding to different signalling networks, all of which are biologically plausible. Parameters in these models were fitted so that they all represented wild type data equally well. The models were then compared to current mutant data and some were invalidated. To discriminate between the remaining models we used ideas from control systems theory to determine efficiently in silico an input profile that would result in the biggest difference in model outputs. However, when we applied this input to the models, we found it to be insufficient for discrimination in silico. Thus, to achieve better discrimination, we determined the best change in initial conditions (total protein concentrations) as well as the best change in the input profile. The designed experiments were then performed on live cells and the resulting data used to invalidate all but one of the remaining candidate models.
Conclusion
We successfully applied our method to chemotaxis in R. sphaeroides and the results from the experiments designed using this methodology allowed us to invalidate all but one of the proposed network models. The methodology we present is general and can be applied to a range of other biological networks.
Background
Understanding the connectivity of signalling pathways within organisms has always been an important challenge in biological research. One approach to address this is to study individual parts in vitro and look at protein localisation, homologies and coexpression in order to elucidate signalling pathway connectivity. Another approach is to use genetics and mutants to attempt to work out the pathway connectivity in vivo.
More recently, systems biology approaches have used quantitative measurements to develop mathematical models that can be used for understanding the properties of biological signalling pathways and their connectivity [1–3]. These models are usually a result of cyclic mathematical model/experiment design iterations, which aim to yield maximum information about the system under study [4, 5]. It is well known, however, that models of comparable complexity corresponding to different pathway connectivities may fit experimental data equally well, leaving the researcher with the dilemma of which model is correct [6]. The question of which model is correct is actually impossible to answer as model validation is a misnomer. The related question of which models are invalid can be answered if appropriate data are available, i.e. the inability of a model to reproduce a data set renders it invalid. Applied in this way, model invalidation can be used to reduce the number of possible models [7], hence narrowing down the number of possible network connectivities.
Previous work in model invalidation emphasised the importance of using timevarying inputs to the system under study. One method of invalidation is to apply a dynamical input and try to maximise the difference in the phase shift of two competing deterministic models [8]. The disadvantage of this approach is that phaseshift is usually difficult to quantify, especially with noisy data. Another approach is presented in [9], where the authors develop a dynamic modelbased controller and an input profile that drives the system output along a prescribed target trajectory. However, this approach requires the implementation of a controller in the laboratory which may prove difficult. Other approaches for model invalidation are presented in [5, 10, 11], and the references therein. These lackoffit methods are used to invalidate models in a statistical manner, but there can be problems with this approach as it relies on large data sets and focuses on obtaining reliable parameter estimates rather than network connectivities. There is also the issue that a wide range of model parameters could give a very similar model output and these methods would have difficulty coping with this.
In this paper we present a new method for developing mathematical models of biological signalling networks, aiming to understand biological network structure. Given a set of experimental data, models corresponding to competing network connectivies are first constructed, all of which can explain wild type data equally well. Then, experiments which maximally discriminate between models corresponding to different networks are designed systematically. These experimental results are used to invalidate models of these networks, resulting in a cyclic process which aims to produce a mathematical model corresponding to a signalling pathway structure which explains all available wild type and mutant experimental data. Our method is applicable to biological pathways for which it is possible to experimentally modify the input profile and measure the output simultaneously. One such system is the chemotaxis signalling pathway, for which tethered cell experiments can be performed in a flow cell and used to measure the response of the system to dynamic ligand concentration profiles.
Chemotaxis is the biasing of movement towards regions of higher concentrations of beneficial or lower concentrations of toxic chemicals by altering the frequency of flagella switching [12]. The signalling pathway within E. coli is well understood and is a simple circuit with one feedback loop [13]. The receptor in the system is a M ethylaccepting C hemotaxis P rotein (MCP) that senses ligands outside the cell. Associated with the MCP is the histidine protein kinase called CheA. Binding of certain ligands to MCP decreases the autophosphorylation rate of CheA. CheA can transfer phosphoryl groups to two possible response regulators, CheY and CheB. CheYP interacts with motor binding sites of the multiple E. coli flagella motors causing a change in direction. The receptors are constantly methylated by the action of a methyltransferase CheR, while CheBP acts as a methylesterase to demethylate the receptor, making it less responsive to ligand binding. This creates a feedback loop, allowing for adaptation. Adaptation allows E. coli to react to changes of the concentration gradient of ligands and not to changes in concentrations per se. Finally, CheZ acts to dephosphorylate CheY and to terminate the signal.
In this paper, we apply our method of model invalidation to different possible and plausible connectivity structures in the chemotaxis pathway of R. sphaeroides. In particular, we use the model invalidation/experiment design cycle outlined above, in order to shed light on the signalling pathway in R. sphaeroides.
Results
A method for discriminating between competing network models
A general method for designing experiments so as to render the outputs of candidate models as different as possible is described. These experiments can then be performed in the laboratory and, when compared to the model predictions, allow the invalidation of some of the candidate models, even when the experimental data set is noisy.
Our method involves the development of ODE models corresponding to different signalling pathway connectivities, all of which can explain current wild type experimental data equally well. The models have in common all currently known interactions and differ in that each model represents a new speculative pathway connectivity. Some parameters in these models are known and some others are unknown. Thus in principle, one can develop a "set" of models with uncertain parameters for a particular signalling structure, each member of this set representing the wild type experimental data equally well. We would want to discriminate between these "sets" of models (which represent signalling pathway structures) but since this is a very hard problem, our method uses a nominal model from each set and then designs an experiment in order to discriminate between these nominal models. Once the discriminatory experiments between the nominal models are designed in silico, we a posteriori assess the discriminatory property of the stimulus by simulating the behaviour for many, randomly chosen, models within these sets in order to see whether the outputs from the two model sets remain distinguishable.
The discriminatory experiments were designed by initially determining the input profile that maximises the magnitude of the squared output difference between models over time. This is typically an optimal control problem, which is often laborious to solve and can result in an input profile that is difficult to realize experimentally [18]. In order to design the input profile in a numerically efficient manner, we used the following result from linear systems theory: for a linear timeinvariant system, the input that produces the largest energy output given an input of fixed energy is a (truncated) sinusoid of a particular frequency [19]. The frequency can be obtained using a Bode plot (see Methods section), a tool which is often used in control and systems theory [20]. Thus, we first linearised the models and then determined this particular frequency of a sinusoidal input corresponding to the error system, i.e. the difference between the two models.
If the above method is insufficient to discriminate between the models then a further method of 'mutating' the two models in biological terms or changing parameters/initial conditions (e.g. alter the initial protein concentrations) in engineering terms is applied to achieve discrimination. These changes can be tested in silico to determine those which will discriminate best between the models under test, before undertaking experiments. The exact nature of the perturbation that can be performed will vary with the system being investigated, but could include altering protein levels by knockout, knockdown (e.g. RNAi in eukaryotes), protein over expression, etc. Thus the space of perturbation which can be searched will be defined by what is possible to implement experimentally in the biological system under investigation.
The designed experiments can then be implemented and the data obtained from such laboratory experiments used to help us differentiate between the models under study, by invalidation. In the above procedure we used deterministic models and a worstcase input design procedure, therefore even if the data are noisy, we expected that we would still be able to invalidate some of these models.
Determining pathway connectivity in R. sphaeroides chemotaxis
We applied our model invalidation method to elucidate unknown interconnections within the R. sphaeroides chemotactic signalling pathway. In vitro phosphorylation measurements have determined some of the internal connectivity and genetic work has shown that only CheY_{6} and either CheY_{3} or CheY_{4} are required for motor switching and that CheY_{6} can bind to the motor [21] (Figure 1). Whilst it has been shown in vitro that all CheYs _{16} can interact with FliM, the rotor switch [22], it is currently unknown whether CheY_{3} and CheY_{4} proteins cause flagella motor switching or whether some CheYs have an effect on other parts of the signalling pathway to influence the motor indirectly, for example, by acting as a phosphate sink.
Model creation and parameter estimation
We obtained the following values for the parameters for the different models (the units for values of K_{1} are in μ M, and in for K_{2} and K_{3}):

blue: K_{1} = 1, K_{2} = 16.5 and K_{3} = 0.015

red: K_{1} = 25, K_{2} = 3250 and K_{3} = 0.075

grey: K_{1} = 1, K_{2} = 97.5 and K_{3} = 0.0023

magenta: K_{1} = 1, K_{2} = 0.4 and K_{3} = 0.011
We assumed that these parameters are the same for the polar and cytoplasic clusters, so .
We used the Pearson productmoment correlation coefficient as a measure of correlation between model prediction and data. For 200 ≤ t ≤ 800 seconds, we found a good correlation between data and the predictions made by the models. For example, we obtained a coefficient of 0.9466 for the blue model and of 0.9544 for the red model.
Reaction rates used in the R. sphaeroides chemotaxis models
reaction  parameter(s)  value(s) 

(R_{1}) A_{2} → A_{2p}  k _{1}  0.03 s^{1} 
(R_{2}) A_{2p+}B_{1} ⇌ A_{2} + B_{1p} 
 0.035(μ Ms)^{1}, 0.01(μ Ms)^{1} 
(R_{3}) A_{2p}+ Y_{3} ⇌ A_{2} + Y_{3p} 
 0.065 (μ Ms)^{1}, 0 
(R_{4}) A_{2p}+ Y_{4} ⇌ A_{2} + Y_{4p} 
 0.004(μ Ms)^{1}, 0 
(R_{5}) A_{2p}+ Y_{6} ⇌ A_{2} + Y_{6p} 
 0.0006(μ Ms)^{1}, 0 
(R_{6}) A_{2p}+ B_{2} ⇌ A_{2} + B_{2p} 
 0.0035(μ Ms)^{1}, 0.01(μ Ms)^{1} 
(R_{7}) B_{1p}→ B_{1}  k _{7}  0 
(R_{8}) Y_{3p}→ Y_{3}  K _{8}  0.08s^{1} 
(R_{9}) Y_{4p}→ Y_{4}  K _{9}  0.02s^{1} 
(R_{10}) Y_{6p}→ Y_{6}  K _{10}  0.1s^{1} 
(R_{11}) B_{2p}→ B_{2}  K _{11}  0.015s^{1} 
(R_{12}) A_{3}A_{4p}+ Y_{6} ⇌ A_{3}A_{4+} Y_{6p} 
 0.1 (μ Ms)^{1}, 0 
(R_{13}) A_{3}A_{4p}+ B_{2} ⇌ A_{3}A_{4} + B_{2p} 
 0.006(μ Ms)^{1}, 0.07(μ Ms)^{1} 
(R_{14})A_{3}A_{4} → A_{3}A_{4p}  k _{14}  0.02s^{1} 
We observed that the difference between the outputs of the two models under consideration was insufficient to allow discrimination between the models experimentally using this discrimination technique. Therefore, we sought experimental perturbations which resulted in a larger difference between the outputs of the different models. Possible experimental perturbations in this system involve deletions of one or more components, over expression or under expression of a protein component and growth of the bacteria in different growth conditions; the latter results in largescale changes of the expression of the chemotaxis operons and hence protein concentrations.
We then used IPTG to experimentally induce the over expression of CheY_{4} from an expression plasmid. The level of CheY_{4} over expression as a population average was verified using semiquantitative western blotting (data not shown). We observed that when CheY_{4} is overexpressed to five times the wild type level the tethered cell trace was similar to that of wild type cells (Figure 6B); this suggested that the model, in which the CheY_{6} and CheY_{3} or CheY_{4} bind to the motor cooperatively (blue model) and do not have other interactions, is invalid, as it is unable to represent our new experimental data (Figure 6B). Since the red model, where there is an interconnection between CheY_{3} and CheY_{4} and the cytoplasmic cluster altering the CheA_{4} kinase activity, can represent this data, it is the only model we cannot invalidate (Figure 6B).
Discussion
We have developed a methodology for elucidating biological signalling pathways which we applied to understanding the chemotactic signalling pathway in R. sphaeroides. Our method differs from many other methods currently in use in that is based on manipulation and observation of the entire biological system under in vivo conditions and, importantly, offers a systematic approach to model invalidation based on cyclic model development/experiment design iterations. In particular, it does not aim at designing experiments for the refinement of parameter values, but rather at identifying possible interconnection structures. By doing so, our method considers a model representing a specific structure that it aims to invalidate.
We applied our method to a real biological system and successfully managed to invalidate some potential signalling pathway network models. Thus, this method helped to more rationally consider the interconnectivity of the chemotactic signalling pathway in R. sphaeroides through relatively simple inputoutput experiments in vivo and helps to design rational future experiments.
The method is applicable to other biological systems but requires an experimental setup where it is possible to control the input and measure the output simultaneously. Even with this limitation the method could be used for both pathway determination and parameter determination, where multiple models with different parameters could be systematically invalidated. Our method could potentially be used to help annotate and understand signalling pathways in nonmodel organisms, using the information from model organisms as a guide for the first model generation step. These could then be tested and models invalidated. A parameter search could also be performed by creating models which deviate from the model organism data and then, through model invalidation, determining which parameter sets are able to fit the experimental data. In all these cases though, the input of the system must be under experimental control and the output easily measured.
A potential limitation would be the case in which two models from different network structures produce outputs that are indistinguishable under all potential experimental conditions. The problem in this case could lie with the particular properties of the network and their discriminatory nature  i.e., whether such differences in structure are identifiable from inputoutput experiments [23]. Performing all calculations in silico before any experiments are undertaken is important in saving experimental time ensuring that experiments are only performed when the results will have discriminating property between the models under test.
In the system studied in this paper the modelling and experiment design was done on average cell traces for R. sphaeroides responses. The reason for this is that single cell behaviour is noisy, and in any case the model parameters we have available at the start of the procedure correspond to the average system behaviour. The total protein concentrations, for example, are the population average determined for these growth conditions. Thus the model output is for the 'average' cell hence we compared the output to the average of the tethered cell data.
When constructing the models we used in vitro rates as it is very difficult, if not impossible to measure these rates in vivo. The multiple homologues in R. sphaeroides would prevent a FRET assay such as the one applied to E. coli being used. Thus the parameters are the best estimate of the real value. To allow for this uncertainty we considered models where these parameters were varied to ensure robust invalidation. Also for simplicity in the modelling we did not consider the increased concentration of the CheAs at specific points in the cell due to clustering nor the effects of clustering. This is because the exact concentration of CheA in specific regions and the effects of clustering are currently impossible to measure experimentally. The fitted parameter K_{1}, which relates the activation of the ligand on the receptor, includes the effect of clustering and as such it is accounted for in our models. To ensure that this is robust we also allowed the concentration of the CheA proteins to vary sufficiently and our invalidation conclusion was still correct.
Interestingly, the connectivity that best fits the experimental data suggests that in R. sphaeroides CheY_{3}P and/or CheY_{4}P do not bind cooperatively with CheY_{6} to FliM. The model we have been unable to invalidate suggests that CheY_{3}P and CheY_{4}P form a link between the polar and cytoplasmic signalling clusters, helping transmit the signal between the two clusters. We have through modelling and comparing to experimental data invalidated a model representing the previously held hypothesis that CheY_{3} and CheY_{4} act only as phosphate sinks for the system. We have also invalidated a model with strict cooperative binding of CheY's at the motor and as such our technique adds to the body of knowledge on R. sphaeroides chemotaxis.
Conclusion
In this paper we have developed a control engineering method for elucidating biological signalling pathways and applied it to a real system. This method is based on multiple model creation and subsequent invalidation using in silico designed experiments. This is a general method that can be applied to other biological pathways where it is possible to control the input and measure the output in simultaneously. We used the method to invalidate all but one model for the chemotaxis signalling pathway in R. sphaeroides and in doing so have invalidated models of certain possible connectivities.
Methods
Strains and growth conditions
Strains and plasmids
Strain/Plasmid  Characteristics  Source 

R. sphaeroides WS8N  Spontaneous nalidixic acid resistant mutant of wild type WS8  [34] 
R. sphaeroides JPA421  WS8N with the cheY_{4} gene deleted by genomic replacement  [35] 
E. coli S171 λ pir  A strain capable of mobilising pAE and pAY4 into R. sphaeroides, Sm^{R}  [36] 
Plasmids  
pIND4  Over expression vector with a lac inducible promoter, capable of replication in R. sphaeroides  [25] 
pIND4Y4  pAE containing cheY_{4}  [29] 
Protein over expression
pIND4Y4 was transformed into S17 and conjugated into R. sphaeroides as described previously [25]. IPTG was added to cell culture at 100 μ M and the cells incubated for ~16 hours at 30°C until OD_{700} is 0.6, where the cells were then used for tethered cell analysis.
Tethered cell analysis
R. sphaeroides was grown to an OD_{700} of 0.6 and then 4 × 1 ml of cells were harvested. 3 × 1 ml were saved for western blot analysis. The remaining 1 ml was pelted and resuspended in tethering buffer (10 mM NaHEPES pH 7.2 containing chloramphenicol at 50 μ g/ml) and incubated at 30°C for 30 mins.
The cells were then tethered in a humidity chamber by incubation of 10 μ l of cell suspension on a coverslip with 1 μ l of antiflagellar antibody for 30 mins.
The coverslip was then loaded onto a flow chamber and tethering buffer passed through for 5 min to remove free cells. After this period the tethered cells were observed under phase contrast at 1000 × magnification. Real time recordings were made on videotape. Tethering buffer with and without Propionate (sodium salt) was passed through the chamber at a rate of 0.09 ml/min.
The video recordings were analysed with the Hobson Bactracker (Hobson Tracking Systems) using the program Arot7. The rotation rate of the cells was measured by detecting the position of the cell every 50 ms. The data obtained were smoothed (100 points), averaged (for as many cells as available  at least 20 per graph) and plotted.
Western blotting
In order to determine protein concentrations semi quantitative western blots were employed as described previously [26].
Modelling the chemotaxis pathway in R. sphaeroides
Our model is split into three modules: sensing, transduction and actuation.
Sensing
 (i)
Receptors can be in different states of methylation. For simplicity, we assume that receptors are either methylated or not.
 (ii)
Only methylated receptors, R^{ m }( ), can be in an active state, R^{ a }( ).
 (iii)
Autophosphorylation of CheA_{2}P (phosphorylation of CheA_{3} by CheA_{4}P) occurs only when the MCP (Tlp) receptor is active.
 (iv)
Binding of the ligand to a receptor inhibits its activity.
 (v)
CheB_{1}P (CheB_{2}P) binds only to active receptors in order to demethylate them.
 (vi)
CheR_{2} (CheR_{3}) binds only to inactive receptors, R^{ i }( ).
 (vii)
The number of CheR_{2} (CheR_{3}) is constant.
We incorporated assumptions (iv), (v) and (vii) into our model as follows:

We let the number of active receptors depend reciprocally on ligand concentrations L ( )  see below for details.

We represented the action of CheB_{1}P (CheB_{2}P) through the following term: K_{2}[R^{ a }][B_{1p}] .

We represented the action of CheR_{2} (CheR_{3}) by the constant reaction rate K_{3}( ).
where ε is a small constant representing disturbances at the input of the cytoplasmic sensing cluster, with nominal amplitude 0.001. We let α = 0.1 and R^{ t }= = μ M.
We assumed that the cytoplasmic sensing cluster senses extracellular ligand concentrations indirectly; for example, could be internalised attractants, a byproduct of the internalisation process or a metabolic response to it. For simplicity, we assumed an affine relationship between L and . Note that for the red model is a function of CheY_{3}P and CheY_{4}P concentration levels, reflecting their effect on the activity of the cytoplasmic sensing cluster.
Transduction
Average copy number of chemotaxis proteins
Number of copies of  value 

CheA_{2}  26000 
CheY_{3}  1000 
CheY_{4}  4000 
CheA_{3}A_{4}  12000 
CheY_{6}  51500 
CheB_{1}  23000 
CheB_{2}  3000 
where [(A_{3}A_{4})_{ i }] denotes the case when CheA_{4}P cannot phosphorylate CheA_{3} due to the action of CheY_{3} and CheY_{4}, and 0.001 ≤ k_{15} ≤ 1; the latter is to say that the behaviour of the magenta model remains virtually unchanged within this parameter range of k_{15}.
We obtained the value of k_{1}, the reaction constant of the autophosphorylation of CheA_{2}, from in vitro experiments in the absence of the influence of receptors. However, when membrane receptors are in their active state they accelerate the autophosphorylation of CheA_{2}. Thus, we modified the reaction constant to ; that is, we assumed that the in vitro reaction rates correspond to the case when receptors are fully active. Similarly, the phosphotransfer from CheA_{4}P to CheA_{3} is accelerated when cytoplasmic receptors are active and we modified the reaction constant to .
The reaction rates were obtained by fitting parameters to data from in vitro experiments [16, 31] (Table 1).
Actuation
We set S = 0.125, which means that saturation occurs at 8 rot/sec.
Parameter fitting
In order to fit model parameters such that the model represents well the experimental data, we minimised the 2norm of the vector, whose entries consist of the errors between data and predictions from the discretised version of the ODE models representing the individual chemical reactions investigated in vitro. The 2norm of vector x is given by , where n is the length of x.
In order to fit any remaining model parameters to data from in vivo experiments, we simulated the model and minimised the 1norm of the vector, whose entries are the errors between data and model predictions [32]. The 1norm of vector x is given by x_{1} + x_{2} + ⋯ + x_{ n }
Fitting model parameters for receptor activation
Parameters K_{1}  K_{3} are unknown and cannot be easily measured by experiments. To determine these we fitted the models to wild type tethering data. We performed tethered cell experiments where we applied 100 μ M of attractant for 5 minutes and then removed it. We then ran a simulation of each model (Figure 3). To obtain K_{ i }s and s, we minimised the error between computer simulations and data following the online fitting procedure. For simplicity and because we are fitting the models to a single model output only, which is contaminated with noise, we let K_{ i }= .
Least squares minimisation to fit phosphotransfer parameters
where p corresponds to the number of measurements.
Using phosphotransfer data [16, 31] and the above method under the constraint that the rates are nonnegative, we obtained values for parameters k_{1}  k_{14} (Table 1).
Online fitting
Because the least squares optimisation is applied to a discretised version of a continuous time model, we applied the method of steepest descent online  that is, we minimise the error between simulations of the continuoustime model and the data  to improve the fit using the values obtained before.
Experiment design
Input design
Here, x_{ i }is the vector with the different states (for example, concentrations of different proteins) of the model i, i = 1, 2, f_{ i }and g_{ i }are functions that represent the different models, and u is the input, which can be externally controlled; for example, ligand concentrations.
Here, τ denotes the duration of the experiment. To obtain an input that maximises (27) for a nonlinear system of the form (24)(25) is difficult. However, if we relax the problem to obtaining the sinusoidal input to the system given by the linearisation of (24)(25) that maximises (26) then it can be solved systematically as we show in the following.
where A, B and C are matrices, whose entries depend on the model parameters, and u(t) = a sin(ωt) is a sinusoidal input with angular frequency ω and fixed amplitude a. System (27) is the so called state space representation of the model in the time domain. It is common in control systems engineering to investigate the behaviour of such a system in dependency of ω, which requires to transform the system to the so called frequency domain. If matrix A is Hurwitz (stable) then after some transient behaviour output y is given by a sinusoidal wave that is, first, phase shifted with respect to u and second, has amplitude . Here we are only interested in the maximum of with respect to ω.
Note that in our analysis the only independent variable is ω. We replaced s by jω (s = jω), where . The Bode magnitude plot shows the value of . Thus, the maximum in the plot provides the frequency that will maximise the output to input ratio.
We linearised (9) to (20) to obtain a system of the form given by (27). We determined the frequency of the sinusoidal input that will result in the largest difference between the model outputs as described above. Finally, we simulated the experimental output of this optimal frequency.
Initial conditions design
Choosing the change or possible combination of changes that provides best discrimination between the different models is a combinatorial problem and difficult to solve efficiently. However, because for our chemotaxis models the range of possible changes was relatively small, due to being limited by what can be implemented experimentally, we performed a brute force search, checking all possibilities. The fivefold over expression of CheY_{4} under microaerobic growth conditions yielded the best result in silico in discriminating between the red and the blue model.
Robustness of invalidation
In order to assess whether our invalidation conclusions are robust to parameter changes, we performed a sensitivity analysis on our models. We ran 500 simulations allowing parameters k_{1} and k_{14} to vary by ± 15% (Figure 4C). Even with these variations, we were able reach the same conclusions.
Declarations
Acknowledgements
This work was funded by the EPSRC, project E05708X, a Royal Society Wolfson Merit Award, and a Praelectorship from Lincoln College, Oxford. The authors would like to thank Dr Marcus Tindall and Dr Steven Porter for useful discussions about this project. The authors would also like to thank the reviewers for their helpful comments.
Authors’ Affiliations
References
 Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387: 913917. 10.1038/43199View ArticlePubMedGoogle Scholar
 Rao CV, Kirby JR, Arkin AP: Design and diversity in bacterial chemotaxis: a comparative study in Escherichia coli and Bacillus subtilis. PLoS Biol. 2004, 2: E49 10.1371/journal.pbio.0020049PubMed CentralView ArticlePubMedGoogle Scholar
 Palsson B: The challenges of in silico biology. Nat Biotech. 2000, 18: 11471150. 10.1038/81125.View ArticleGoogle Scholar
 Clemens K, Jens T: Systems biology: experimental design. FEBS Journal. 2009, 276: 923942. 10.1111/j.17424658.2008.06843.xView ArticleGoogle Scholar
 Maksat A, Yves FN, Jaap AK, Joke GB: Systems biology: parameter estimation for biochemical models. FEBS Journal. 2009, 276: 886902. 10.1111/j.17424658.2008.06844.xView ArticleGoogle Scholar
 Alon U: An introduction to systems biology: design principles of biological circuits. 2006, Boca Raton, Fla; London: Chapman & Hall/CRCGoogle Scholar
 Anderson J, Papachristodoulou A: On validation and invalidation of biological models. BMC Bioinformatics. 2009, 10: 132 10.1186/1471210510132PubMed CentralView ArticlePubMedGoogle Scholar
 Kremling A, Fischer S, Gadkar K, Doyle FJ, Sauter T, Bullinger E, Allgöwer F, Gilles ED: A benchmark for methods in reverse engineering and model discrimination: Problem formulation and solutions. Genome Res. 2004, 14: 17731785. 10.1101/gr.1226004PubMed CentralView ArticlePubMedGoogle Scholar
 Apgar JF, Toettcher JE, Endy D, White FM, Tidor B: Stimulus design for model selection and validation in cell signaling. PLoS Comput Biol. 2008, 4: e30 10.1371/journal.pcbi.0040030PubMed CentralView ArticlePubMedGoogle Scholar
 Chen BH, Asprey SP: On the design of optimally informative dynamic experiments for model discrimination in multiresponse nonlinear situations. Ind Eng Chem Res. 2003, 42: 13791390. 10.1021/ie0203025.View ArticleGoogle Scholar
 Faller D, Klingmuller U, Timmer J: Simulation Methods for Optimal Experimental Design in Systems Biology. SIMULATION. 2003, 79: 717725. 10.1177/0037549703040937.View ArticleGoogle Scholar
 Baker M, Wolanin P, Stock JB: Signal transduction in bacterial chemotaxis. Bioessays. 2006, 28: 922. 10.1002/bies.20343View ArticlePubMedGoogle Scholar
 Wadhams GH, Armitage JP: Making sense of it all: Bacterial chemotaxis. Nat Rev Mol Cell Biol. 2004, 5: 10241037. 10.1038/nrm1524View ArticlePubMedGoogle Scholar
 Sourjik V, Schmitt R: Phosphotransfer between CheA, CheY1, and CheY2 in the chemotaxis signal transduction chain of Rhizobium meliloti. Biochemistry. 1998, 37: 23272335. 10.1021/bi972330aView ArticlePubMedGoogle Scholar
 Wadhams GH, Warren AV, Martin AC, Armitage JP: Targeting of two signal transduction pathways to different regions of the bacterial cell. Mol Microbiol. 2003, 50: 763770. 10.1046/j.13652958.2003.03716.xView ArticlePubMedGoogle Scholar
 Porter SL, Armitage JP: Phosphotransfer in Rhodobacter sphaeroides Chemotaxis. J Mol Biol. 2002, 324: 3545. 10.1016/S00222836(02)010318View ArticlePubMedGoogle Scholar
 Porter SL, Wadhams GH, Armitage JP: Rhodobacter sphaeroides: complexity in chemotactic signalling. Trends Microbiol. 2008, 16: 251260. 10.1016/j.tim.2008.02.006View ArticlePubMedGoogle Scholar
 Bryson AE, Ho YC: Applied optimal control: optimization, estimation, and control. Revised printing edn. 1975, New York; London: HemisphereGoogle Scholar
 Doyle JC, Francis BA, Tannenbaum A: Feedback control theory. 1992, New York Toronto: Macmillan Pub. Co: Collier Macmillan Canada; Maxwell Macmillan InternationalGoogle Scholar
 Franklin GF, Powell JD, EmamiNaeini A: Feedback control of dynamic systems. 2005, Upper Saddle River, NJ: Pearson Prentice Hall, 5Google Scholar
 Porter SL, Wadhams GH, Martin AC, Byles ED, Lancaster DE, Armitage JP: The CheYs of Rhodobacter sphaeroides. J Biol Chem. 2006, 281: 3269432704. 10.1074/jbc.M606016200View ArticlePubMedGoogle Scholar
 Ferre A, de la Mora J, Ballado T, Camarena L, Dreyfus G: Biochemical study of multiple CheY response regulators of the chemotactic pathway of Rhodobacter sphaeroides. J Bacteriol. 2004, 186: 51725177. 10.1128/JB.186.15.51725177.2004PubMed CentralView ArticlePubMedGoogle Scholar
 August E, Papachristodoulou A: A new computational tool for establishing model parameter identifiability. J Comput Biol. 2009, 16: 875885. 10.1089/cmb.2008.0211View ArticlePubMedGoogle Scholar
 Sistrom WR: A requirement for sodium in the growth of Rhodopseudomonas spheroides. J Gen Microbiol. 1960, 22: 778785.View ArticlePubMedGoogle Scholar
 Ind AC, Porter SL, Brown MT, Byles ED, de Beyer JA, Godfrey SA, Armitage JP: InducibleExpression Plasmid for Rhodobacter sphaeroides and Paracoccus denitrificans. Appl Environ Microbiol. 2009, 75: 66136615. 10.1128/AEM.0158709PubMed CentralView ArticlePubMedGoogle Scholar
 Martin AC, Nair U, Armitage JP, Maddock JR: Polar Localization of CheA2 in Rhodobacter sphaeroides Requires Specific Che Homologs. J Bacteriol. 2003, 185: 46674671. 10.1128/JB.185.16.46674671.2003PubMed CentralView ArticlePubMedGoogle Scholar
 Alon U, Surette MG, Barkai N, Leibler S: Robustness in bacterial chemotaxis. Nature. 1999, 397: 168171. 10.1038/16483View ArticlePubMedGoogle Scholar
 Gould M: Chemotaxis gene expression in Rhodobacter sphaeroides WS8N. DPhil Thesis. 2006, University of OxfordGoogle Scholar
 Brown M: Control of the Unidirectional Motor in Rhodobacter sphaeroides. D Phil Thesis. 2009, University of OxfordGoogle Scholar
 Armitage JP, Evans MCW: Control of the protonmotive force in Rhodopseudomonas sphaeroides in the light and dark and its effect on the initiation of flagellar rotation. Biochim Biophys Acta. 1985, 806: 4255. 10.1016/00052728(85)900805.View ArticleGoogle Scholar
 Porter SL, Armitage JP: Chemotaxis in Rhodobacter sphaeroides requires an atypical histidine protein kinase. J Biol Chem. 2004, 279: 5457354580. 10.1074/jbc.M408855200View ArticlePubMedGoogle Scholar
 August E, Papachristodoulou A: Efficient, sparse biological network determination. BMC Syst Biol. 2009, 3: 25 10.1186/17520509325PubMed CentralView ArticlePubMedGoogle Scholar
 Zhou K, Doyle JC, Glover K: Robust and optimal control. 1996, PrenticeHall, IncGoogle Scholar
 Sockett RE, Foster JCA, Armitage JP: Molecular biology of the Rhodobacter sphaeroides flagellum. FEMS Symp. 1990, 53: 473479.Google Scholar
 Shah DS, Porter SL, Martin AC, Hamblin PA, Armitage JP: Fine tuning bacterial chemotaxis: analysis of Rhodobacter sphaeroides behaviour under aerobic and anaerobic conditions by mutation of the major chemotaxis operons and cheY genes. EMBO J. 2000, 19: 46014613. 10.1093/emboj/19.17.4601PubMed CentralView ArticlePubMedGoogle Scholar
 Penfold RJ, Pemberton JM: An improved suicide vector for construction of chromosomal insertion mutations in bacteria. Gene. 1992, 118: 145146. 10.1016/03781119(92)90263OView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.