Discussion on the models

Accuracy

We address three major issues regarding accuracy of these models.

Numerical accuracy

The exponential Euler formulation was used for integration. This is an explicit method based on the Forward Euler method, but it assumes exponential decays of variables. This makes it well suited for handling kinetic simulations. The accuracy improves with shorter simulation time-steps.
Implicit simulation techniques with variable time-steps were also examined, but these did not offer significant improvement in simulation times. The simpler Exponential Euler method was therefore used in all the model calculations.

The numerical accuracy of the computation was verified in three ways:

1. Simple kinetic schemes which could be calculated analytically were modeled and simulated results compared to the analytical results.
2. The law of mass conservation and microscopic reversibility principles were used as tests of accuracy in complex reaction schemes.
3. The same model was run at different time steps and the resulting simulated values were compared. The calculations were shown to be stable and to provide convergent solutions for the range of time steps used in the study.

In each case the simulations met accuracy criteria of 5% or better. We did not seek to improve accuracy further, as the experimental uncertainty was considerably larger.

Computations were carried out using the well stirred assumption, i.e. that all molecules had equal access to all others. This means that diffusion and compartmentalization were neglected (see discussion below).

Experimental accuracy and limited data

In carrying out simulations of natural processes, there are two major issues that must be addressed with regards to the experiments. The first is whether all the relevant experimental variables have been considered. The second is whether the experimental data is accurate enough.
In our models, missing experimental variables would correspond to signaling interactions which we have omitted. There are certainly many such missing interactions. This is partly due to our inability to model every reported interaction, and partly because in all likelihood further interactions remain to be discovered. We would, however argue, that we have probably accounted for many of the important interactions. Current experimental strategy is well suited to defining major interactions first, and then homing in on details. For instance, the proliferation of isoforms of key enzymes is a detail which we have represented by including "averaged" models for enzymes such as PKC.

The second issue - experimental accuracy - is one which we have taken pains to estimate (Figure 2F in the paper). We feel that barring the conceptual limitations mentioned below, our simulations are likely to behave in qualitatively and quantitatively very similar ways within wide ranges of experimental uncertainty.

One of the important functions of a theoretical analysis of experimental results is to indicate where more rigourous experimental data would be required to fully understand a system. Biochemical reports are often unusually good from a modellers perspective (Bower JM and Koch, C, Trends Neur. Sci. 15:11 1992.) in having a well-defined set of parameters such as Vmax and Km which are reported complete with estimates of error. However, a few recurring gaps frustrate accurate modelling, even though the raw data should be available from the same experiments. The most prominent is lack of information about time-courses, when reporting Kd for an equilibrium reaction. Direct time-course values such as profiles of reaction progress were available in most cases. In some cases the time-course had to be inferred from reports of the time of onset of a response, or by time-courses of down-stream effects. Additional data to estimate the ratio k2:k3 for Michaelis-Menten enzymes would also be very useful in constraining the models. In our simulations the ratio 4 was chosen as it results in a low percentage of enzyme sequestered in a complex with the substrate over most of the physiological range of enzyme and substrate concentrations. We varied this ratio over an order of magnitude in either direction and observed very little change in simulation results.

Conceptual limitations

Conceptual limitations to our models include the following issues which are topics of active research in many laboratories:

1. Cytoskeletal interactions and reaction channeling. Protein scaffolds have recently been identified as important determinants of reaction specificity in several yeast MAP Kinase cascades. Such cytoskeletal anchors for active molecules can result in an assembly-line organization of signaling reactions that may be quite distinct from solution interactions both in terms of specificity and rates. We anticipate that such interactions may significantly change our understanding of signaling processes, and will need considerably more experimental data and new modeling approaches.
2. Single-molecule interactions. At physiological concentrations, there may be only a few individual molecules interacting in a dendritic spine or at a vesicle release site. In such situations the common approximation of molecular concentration as a continuous variable breaks down, and stochastic Markov processes are required to describe chemical interactions.

Robustness

The modelling approach outlined here is obviously limited by the reliability of experimental information. Experimental error and variability between preparations introduces uncertainty in parameter estimates. Nevertheless, many aspects of our models and probably signalling pathways in general, are quite insensitive to the exact parameter values (Figure 2F). The concentration of Protein Kinase C was one of the most sensitive parameters, allowing 25% variation above or below the reference value while retaining the property of bistability. Most other parameters were robust within a factor of 2 in either direction. Work in preparation indicates that the bistable property of the system exists over a remarkably large region of parameter space. We cannot prove bistability with the current data, which has a wide range of experimental values and large errors. Nevertheless, we feel that the available data suggest that the system exists in a region where bistability is likely.

Diffusion and compartmentalization

All reactions were modeled as point reactions, using the "well-stirred cell" assumption. Spatial aspects of calcium fluxes were examined in a limited manner, using explicit compartmentalization in the calcium model. We represented 3 Ca compartments: extracellular, intracellular and intracellular stores. We also carried out simple 1-dimensional Ca diffusion calculations in the dendritic spine model. Several of our enzymes were also represented as transiting between cytosolic and membrane-bound forms. In our models this did not involve explicit spatial compartments. Instead we limited the reactions which the enzymes could participate in. This treatment clearly oversimplifies the enormous implications of spatial segregation and messenger diffusion in real cells. Nevertheless, a point model such as ours is an essential first step towards building a more comprehensive model which incorporates diffusive and cytoskeletal effects.