Stochastic and deterministic multiscale models for systems biology: an auxin-transport case study

Background Stochastic and asymptotic methods are powerful tools in developing multiscale systems biology models; however, little has been done in this context to compare the efficacy of these methods. The majority of current systems biology modelling research, including that of auxin transport, uses numerical simulations to study the behaviour of large systems of deterministic ordinary differential equations, with little consideration of alternative modelling frameworks. Results In this case study, we solve an auxin-transport model using analytical methods, deterministic numerical simulations and stochastic numerical simulations. Although the three approaches in general predict the same behaviour, the approaches provide different information that we use to gain distinct insights into the modelled biological system. We show in particular that the analytical approach readily provides straightforward mathematical expressions for the concentrations and transport speeds, while the stochastic simulations naturally provide information on the variability of the system. Conclusions Our study provides a constructive comparison which highlights the advantages and disadvantages of each of the considered modelling approaches. This will prove helpful to researchers when weighing up which modelling approach to select. In addition, the paper goes some way to bridging the gap between these approaches, which in the future we hope will lead to integrative hybrid models.


Background
Biological systems are naturally multiscale and to understand their behaviour fully we must understand the interaction of a number of processes that may occur on diverse temporal and spatial scales. To gain insight into such multiprocess and multiscale systems, there is a range of modelling frameworks that could potentially be employed. Different modelling approaches serve to highlight certain aspects of a biological system, and which modelling approach is most appropriate depends on the biological questions that are being addressed, as well as on the available data that could be used to calibrate or validate a given model. In this paper, we present several modelling approaches and show how these can be used to gain understanding of a realistic multiscale systems biology problem. We compare the different modelling approaches to each other and discuss their applicability.
To compare the modelling approaches, we focus on a particular case study: the transport of the hormone auxin through a file of plant cells. Auxin plays a major role in many aspects of plant growth and development [1]. It moves through the plant in a polar manner due to non-uniform spatial distributions of active influx and efflux carriers on the cell membranes [2], and the resulting auxin distributions influence a wide range of processes, including organ initiation [3][4][5][6], vein formation [7][8][9][10][11][12] and gravitropism [13]. Modelling auxin transport is thus an active research area in plant systems biology. The models are inherently multiscale, as cell-scale processes lead to tissue-scale phenomena. To date, the majority of modelling in this area computes solutions by simulating large systems of deterministic ordinary differential equations [2][3][4][5][6][7][8][9][11][12][13][14][15][16][17], and there are relatively few examples of alternative modelling techniques [18][19][20].
This paper complements such previous work by highlighting the benefits of using multiple modelling techniques to gain a more comprehensive understanding of a biological system, in this case auxin transport in plant tissue.
Mathematical modelling is routinely used to study biological phenomena quantitatively, often by describing the dominant physical processes using systems of coupled differential equations and solving these governing equations using analytical and numerical methods; such techniques have been used to study a diverse range of biological processes, including population dynamics, pattern formation, neuron firing and physiological flows (see [21] and [22] and references therein). Mathematical models have denotational semantics in that they represent relationships between quantities as systems of equations. In contrast, computational, or executable, models have operational semantics, and define rules that describe how the modelled system moves from one state to the next [23]. Computational models are executed in the sense that, starting from an initial state, a procedure (in our case a stochastic simulation algorithm) determines the next reaction to apply. This reaction is then applied, giving a new system state, which is then used by the simulation algorithm to determine the next rule to apply, in an iterative procedure. Stochastic processes, unlike their deterministic counterparts, involve an indeterminancy in the evolution of the state of the system. For large numbers of molecules, this stochasticity may be averaged out, giving what appears to be a deterministic process; however, when a small number of molecules is involved, stochastic effects become evident, and in such cases the system may behave in a markedly different way. The inherent noise present in all biological systems is explicitly modelled in discrete stochastic models and can have profound effects on system dynamics, producing behaviour, even for large numbers of molecules, which is markedly different from that predicted by continuous deterministic models; see, for example, [24].
After summarising the biological abstraction which forms a common basis for the models presented in this paper, we describe a stochastic computational model based around a P system framework. Here the number of auxin molecules in each compartment evolves according to rules that move molecules from one compartment to the next. We then describe a deterministic mathematical model in which the auxin concentrations are described by a system of coupled ordinary differential equations. In the deterministic case, we produce two solutions: i) analytical solutions, derived using multiscale asymptotic approaches, and ii) numerical solutions, as is typical in current auxin-transport modelling.

Biological Abstraction
To investigate the benefits of deterministic and stochastic modelling approaches, we focus on a model of auxin transport. Specifically, we model a standard experiment that is used to determine auxin velocities: radio-labelled auxin is added to a source agar block at one end of a segment of stem tissue, auxin then travels through the stem segment, and experimentalists measure the amount of auxin collected in a final agar block at the other end of the stem segment (see [25][26][27] and references therein).
As shown in Figure 1, we model the stem segment as a single two-dimensional line of N cells, where each cell contains a cytoplasm (with length l and width w), and there is a layer of apoplast between neighbouring cytoplasms (which has thickness l and width w). Hence the total length of the stem segment is L = Nl + (N + 1)l. The model assumes that the auxin concentration within each compartment is uniform, which is a reasonable assumption given the small size of the compartments and the rapidity of auxin diffusion.
We suppose that the two agar blocks are rectangular, with length L s and width w. We denote the number of auxin molecules in the source agar block by S n (t) and Figure 1 Model of a single file of cells. In the model, auxin molecules are initially in the source agar block, travel through the stem segment and are collected in the final agar block. We model the stem segment as a single two-dimensional line of N cells for simplicity, and suppose that efflux carriers are located on the downstream face of each cell membrane. We solve for the auxin concentrations in the two agar blocks, the N cytoplasm compartments and the N + 1 apoplast compartments. the number in the collecting agar block by F n (t), where the superscript n emphasises that these quantities are numbers of molecules; the concentrations (i.e. number per unit area) in the two agar blocks are then given by S(t) = S n (t)/(L s w) and F(t) = F n (t)/(L s w), respectively. We suppose that auxin diffuses within the agar blocks with diffusion coefficient D. Auxin in the source agar block diffuses into the adjacent root apoplast region; it then travels through the line of N cells and diffuses from the final apoplast region into the collecting agar block. We denote the number of molecules in the cytoplasm by  Figure 1). Auxin exists in planta in either an anionic or a protonated form. Following previous auxin-transport models, we assume that auxin molecules rapidly associate or dissociate, so that the proportion of these two forms are in equilibrium and are determined by the pH of the region in which it is located and by the auxin dissociation constant, pK. Using the subscripts c and a to refer to the cytoplasm and apoplast respectively, the fractions of protonated and anionic auxin are given by [14] Protonated A pH j pK j     1 In the line of cells, auxin moves between the apoplast and cytoplasms by crossing the cell membrane. The flux of protonated auxin across the membrane, J diff , is passive, whereas the flux of anionic auxin is mediated by PIN efflux carriers that are present on the downstream face of each cell membrane (see Figure 1); following [13], we model the anionic flux across the cell membrane, J PIN , using Goldman-Hodgkin-Katz theory (see [21] for details). Thus, the components of the flux from the cytoplasm c i to the apoplast a i are given by where and P diff is the membrane permeability of protonated auxin, P PIN is the membrane permeability of anionic auxin, and the dimensionless constant j ≡ -F D V/RT where F D is the Faraday constant, V is the membrane potential, R is the gas constant, and T is the temperature.
At time t = 0, all the auxin molecules are in the source agar block; we prescribe S n (0) = C, and let c i n ( ) 0 = a j n ( ) 0 = F n (0) = 0 for i = 1, 2, ... N and j = 0, 1, 2, ... N. We model a closed system, that is, we assume that, during the subsequent auxin movement, no auxin molecules enter or leave the system. Table 1 summarises estimates of the model parameters based on values reported in the biological literature; these are discussed further in the "Biological Parameter Estimates" section of Additional file 1.

Stochastic Computational Model
To obtain stochastic solutions, we use a multi-compartment stochastic P system framework [28]. Individual molecules of auxin are modelled as objects which move between compartments according to a set of rules associated with each compartment. Compartments of the same type have the same set of rules associated with them, and each rule has an associated stochastic reaction constant k which determines the rate at which the rule transports molecules (see Table  2). We define rules for auxin diffusion between source/collecting blocks and the apoplast, membrane diffusion from apoplast to cytoplasm, and carriermediated efflux from cytoplasm to apoplast. These reaction constants are related to the parameters given in Table 1 via diffusion form agar blocks to apoplast (further information and derivations can be found in the "Model Derivation" section of Additional file 1). We note that in this framework we set the small parameter A c  to zero to enable efficient execution. The computational model is executed using a novel multi-compartment Monte Carlo stochastic algorithm using the mcss simulator, which is part of the Infobiotics workbench, a freely available software suite for designing, simulating and analysing multiscale executable systems and synthetic biology models [29]. Stochastic algorithms and software supporting multiple compartments have been developed by several research groups [30][31][32][33]. A key difference between these algorithms and our software suite is that our tools support the rapid prototyping of models by facilitating the abstraction of commonly occurring motifs (e.g. regulatory or signalling motifs [34]) with model templates and modules. This, coupled with the facility to explicitly specify a tissue geometry, permits the seamless exploration of "what if " scenarios during model building. For example, to reproduce the behaviour of a mutant which does not produce a particular protein, all that needs to be done is to set the rate constant of the reaction producing the protein to zero. Or, for example, to remove a particular regulatory mechanism, we can simply remove the module representing this regulatory mechanism from the list of modules employed by a particular compartment. The reader is referred to [28] for an in-depth description of the stochastic simulation algorithm we employ, which we briefly summarise here. Essentially, the simulation algorithm determines, within each of the simulated compartments, which rule to apply next and when to apply it. This highlights the key difference between mathematical and computational (or algorithmic) models -mathematical equations describe the change in the values of variable as a system moves from one state to the next, while computational models expose how and why this state change occurs [23,35]. Each run of the simulation gives one possible trajectory of the model through state space. Hence, as well as using individual runs to examine the reasons for state changes, we can execute the model a number of times to estimate the average system behaviour, analyse the system's variability, and identify potential extreme behaviours.
Additionally, the computational model enables model checking, the formal verification of model behaviour, which allows the identification of general biological principles which underlie the observed behaviour of the model [36][37][38][39]. Using algorithms such as those presented in [40], the computational model also facilitates model parameter and structural optimisation, allowing incomplete biological information regarding, for example, model parameter settings and structure, to be established automatically.

Deterministic Mathematical Model
In the deterministic framework, the auxin concentration in each region is governed by the following system of equations Molecules A of auxin are transported from compartment i to a neighbouring compartment at rates k which depend on the type of transport process. For apoplasts, rules R 4 and R 5 are assigned to apoplasts bordered on both sides by cytoplasms. Rule R 7 replaces rule R 5 for the first apoplast, and R 8 replaces R 4 for the final apoplast, which are bordered by the source and collecting agar block respectively.
We now use asymptotic methods to derive an analytical solution for these governing equations, (6)(7)(8)(9)(10)(11): these explicitly contain at least two distinct length scales, namely the cell length scale, l, and the (much larger) tissue length scale, L = Nl + (N + 1)l, -it is this multiscale aspect that we exploit in deriving a continuum formulation. For comparison, we also solve the deterministic governing equations numerically using Matlab's ordinary differential equation solver, ode45. To obtain an analytical solution, we consider the dynamics on the time scale of active transport the length of the tissue, L/P PIN , and nondimensionalise using ( , , , ) ( , , , ), .
The model then depends on the dimensionless parameters To make analytical progress, we construct an asymptotic solution in the biologically appropriate limit in which is small, i.e. the length scale L of the tissue is much larger than that of a single cell, l + l, expanding the concentrations as standard perturbation series. Based on the biologically relevant parameter estimates in Table 1, we assume that the relative length of the agar blocks,L s = O(1) as 0 + , and rescale the small parameterŝ ,ˆ, , where  , , , D A B c a  and P diff are O(1) as 0 + . In the stem segment, the result is a continuum limit. Letting x measure the length along the tissue, such that a(x, t), where x = 0 corresponds to the upper face of cytoplasm i = 1.
Based on these assumptions, we can derive formulae for the leading-order auxin concentration in each compartment (see the "Derivation of Asymptotic Solution" section of Additional file 1 for details). We identify two time scales. On the transport time scale,t = O(1), the source-agar-block concentration does not deplete significantly,Ŝ ≈ 1, and auxin travels through the stem segment with a defined front according to a wave equation, The sink-agar-block concentration is small, and given byˆ(ˆ)ˆˆ.

F t DA a P diff
On a longer time scale,  t t O    ( ) 1 , the cytoplasm concentrations in the stem segment are uniform, and the agar-block concentrations are given by therefore, diffusion within the agar blocks determines the rates at which the source-agar-block concentration depletes and the sink-agar-block concentration increases.

Auxin Concentrations
For both the deterministic and stochastic frameworks, the relative auxin concentrations in the source agar block, collecting agar block, cytoplasm and apoplast regions are shown in Figure 2. For the stochastic model, we show both the mean concentrations and the 95% confidence intervals (which are calculated using 10,000 runs). To compare quantitatively the solutions of the deterministic and stochastic models, we compute the time taken for the agar-block concentrations to reach half their steady-state levels. These results, given in Table 3 and Figure 2, demonstrate excellent agreement between the three modelling approaches.
Computational models, because of their algorithmic specification, are amenable to model checking. Model checking allows formal verification that a model satisfies a prescribed property. Properties are propositions about the state of the model, for example, the amount of species A that reaches a level x. For the computational model, we used the Infobiotics workbench to perform probabilistic model checking and formally determine the probability of the agar-block concentrations reaching half their steady-state levels. Due to computational constraints, we check a reduced version of the computational model (see the "Model Checking" section of Additional file 1 for details). Figure 3 shows that, with 95% confidence, after 215 minutes both the agar-block concentrations will have reached half their steady state, and that the source-agar-block will reach its steady-state concentration slightly before the collecting agar block. To characterise when the deterministic and stochastic models agree and differ, we Concentrations are given relative to the initial source concentration C/(L s w). For the stochastic solutions, the mean concentrations are calculated over 10,000 runs, and 95% confidence intervals are given as grey ranges. Time taken to reach half the steady-state concentrations for the deterministic and stochastic models. 95% confidence intervals are also given for the stochastic model. subtract the concentration predicted by the analytical deterministic model from the mean concentration predicted by the stochastic model and divide this value by the standard error given by the stochastic model; the predicted concentrations are then considered to differ significantly if the absolute value of the result is greater than two, that is, if the difference between the two predicted concentrations is greater than two standard errors. Figure 4 shows that there is no significant difference between the source concentrations predicted by the deterministic and stochastic models. In contrast, in the collecting agar block there is a significant difference for the initial 196 seconds, with a maximum difference of six standard errors. The deterministic approach takes the concentrations to be continuous while the stochastic model considers individual molecules. As a result, for the first 138 seconds, the stochastic model typically predicts a zero concentration, whereas the deterministic model predicts a small positive concentration. As shown in Figure 4, for the remaining time, although not significantly different, the stochastic model consistently predicts larger concentrations than the deterministic one.
To assess how the initial source concentration affects the variability of the auxin concentrations, we performed 10,000 runs of the stochastic computational model with initial concentrations of 0.1 nM and 1 nM (see Table 4). The mean variability and its standard deviation in all compartment types decreases by around 69% as the initial concentration is increased from 0.1 nM to 1 nM. Theoretically, the amount of noise in a stochastic simulation is of the order of n , where n is the number of molecules [41]. Since 0 1 1 . / ≈ 0.31, we would theoretically expect a 69% decrease in noise, in agreement with the simulation results.

Model applicability
Selecting the appropriate modelling approach for a given problem involves a number of factors including a researcher's judgement of the time available for model development and the computational resources available. A key factor is the understanding they wish to obtain, for example, arriving at the best possible model for a given system (in which case model development time might not be too important), or generating alternative plausible models representing competing experimental hypothesis (in which case the availability of rapid prototyping, averages and outlier behaviours are needed). Although difficult to give all-encompassing guidelines, in this section we briefly discuss general practical considerations arising from our case study that are involved in determining an appropriate modelling approach. The stochastic computational model needs to be executed many times to assess the average behaviour of the system for a given set of parameter values, and therefore the computational cost of solving the stochastic model will generally be greater than that of the deterministic model. The execution time of the computational model depends on the values assigned to the reaction constants and the initial numbers of molecules present in the system (as with more molecules there are more reactions per second). Table 5 summarises the time taken to execute the computational model for several initial concentrations, and compares these execution times with the time taken to solve the deterministic numerical model. Tests were performed on a single processor (AMD Athlon 64 X2 Dual Core 5600+ 2.8 GHz with 1 GB of memory) -using multiple processors would reduce these execution times. As shown in Table 5, as expected, the execution time of the stochastic model increases as the initial number of molecules increases. For initial auxin concentration above 1 nM, the stochastic model will take a considerable amount of time, and therefore in this case it would be advisable to use a deterministic approach. In contrast, the asymptotic solution of the mathematical model does not involve any computational cost (other than the negligible time required to evaluate the derived expressions at a given timepoint). However, the computational execution time does not reflect the true time taken to obtain the solution, as deriving the asymptotic formulae may take longer than producing the numerical or stochastic models.
When considering the time costs of different modelling approaches, a key consideration is the number of different parameter choices one wishes to investigate. The formulae from the analytical approach clearly show how the parameter estimates affect the predicted concentrations and transport speeds provided the scaling, (14), holds. However, we would need to derive a new asymptotic solution if we wanted to consider different parameter regimes. In addition, the asymptotic method presented here is can only be applied if we can consider the tissue to be a continuum, which is only appropriate if the rate of transport between the cells is not too small [21]. Using the numerical deterministic and stochastic approaches, one can use any parameter values in the simulations by making simple changes to the numerical code without changes to the underlying model. However, using these methods, we would need to execute many simulations to thoroughly understand how the dynamics are affected by the parameter values, which needs to be balanced with the increased execution time for the stochastic model for certain parameter regimes.
In summary, the asymptotic model is applicable only to the specific parameter regime for which it was derived, but allows rapid evaluation of the behaviour of the model within these bounds. However, to explore model behaviour outside of the given parameter regime, the asymptotic solution will need to be derived anew. The stochastic computational model allows any parameter regime or spatial scaling to be explored without further reformulation of the model and formally captures the mechanisms involved in producing a given phenotype. For some regimes the execution time of the model will be considerable, although this time can be ameliorated through the application of more computational resources or parallel computation.

Auxin-Transport Speed
Auxin-transport experiments aim to investigate the movement of auxin through plant tissue. We have modelled an experimental protocol, as described in [25,27], that has been used to consider both the distance moved by auxin molecules per unit time (the velocity) and the amount of auxin passing through the tissue per unit time (the flux). In the deterministic model, auxin from the source agar block moves through the tissue with a defined front, and the asymptotic solutions provide a simple formula for the speed of transport (15). However, in practice there will be stochasticity in the auxin  movement. The stochastic model predicts when the first auxin molecule appears in the collecting agar block, and the transport speed can be calculated by dividing the total length of the stem segment, L, by the average time taken for the first molecule of auxin to appear in the collecting agar block (this time is calculated by averaging over 10,000 runs).
The results, summarised in Table 6, show that the stochastic approach predicts a greater auxin-transport speed (3.38 ± 0.3 cm·h -1 ) than the asymptotic deterministic solution (1.95 cm·h -1 ). We note that this discrepancy is in agreement with the discussions in "Auxin Concentrations" section above where we showed that the different modelling frameworks predicted significantly different concentrations in the collecting agar block at early times.
Auxin velocities are generally thought to be around 1 cm·h -1 , which is fairly close to our predictions, given that the velocity depends on the parameter estimate for P PIN and this value is not well characterised. One reason for the difference between the experimental auxin velocity and the model predictions may be differences in auxin detection sensitivity between the wet experiments and models. The stochastic model enables us to predict the time at which the first molecule of auxin enters the collecting agar block. However, in the wet experiments, a certain amount of auxin must accrue in this agar block before detection is possible. We can estimate the amount of auxin present in the collecting block from our models. If we consider an experiment with 12, 044 molecules of auxin and assume that, in line with experimental results, the auxin-transport speed is 1 cm·h -1 , then the time taken for auxin to travel the length of the stem segment is 723.78 s (0.20105 hr). The deterministic model gives the number of molecules in the collecting agar block at this time to be 191 and 246 from the numerical and asymptotic solutions respectively (we note that this accuracy is within the expected range for the asymptotic solution). The mean number of auxin molecules at this time in the collecting agar block calculated over 10,000 runs of the stochastic model is 299.40 ± 15.84. Thus, to determine accurately the presence of auxin in the collecting agar block, the experimental apparatus used must have a sensitivity of 1.6 pM for an agar block of the same size as we simulated, and a finer resolution for larger agar blocks.
As discussed in [27], the majority of auxin-transport measurements report the flux of auxin transport rather than the auxin velocity and so consider the amount of auxin that has moved through a specified distance of tissue in a constant amount of time. However, the asymptotic solutions of the deterministic model demonstrate that diffusion within the agar blocks may significantly affect the auxin concentration within the collecting agar block, and therefore the auxin fluxes measured. The analysis presented in the "Derivation of Asymptotic Solution" section of Additional file 1 shows that there are two disparate time scales: on a short time scale, auxin is transported through the stem segment, whereas over longer ones, the auxin concentrations are almost uniform throughout the stem segment, and the dynamics are dominated by diffusion within the agar. It is clearly important to be aware of these two processes when interpreting experimental results. If an auxintransport experiment were carried out over several hours, the auxin concentration in the collecting agar block would be determined predominantly by the diffusion rate. We emphasise that these conclusions are based on the assumption that the agar-block length is comparable to the stem-segment length -the effect of diffusion within the agar blocks will be less significant with smaller agar blocks.

Conclusions
In systems biology, models are typically deterministic and a biological problem is translated into large systems of ordinary differential equations that are solved numerically. However, this is not the only option, and in this paper we have demonstrated three different modelling approaches: (i) deterministic numerical; (ii) deterministic asymptotic; and (iii) stochastic computational. As expected, particularly given that the dynamics can be described by a system of linear governing equations, there is excellent agreement between the three methods.
We have focussed our case study on auxin transport, as this is inherently multiscale with cell-scale dynamics creating the tissue-scale phenomena of interest. The numerical, deterministic method focusses on computing the cell-scale dynamics, whereas the asymptotic method makes use of the multiscale nature of the system: in the asymptotic results, we consider the auxin concentrations on the cell scale, and exploit the relatively small dimensions of the cells to determine how the cell-scale dynamics lead to effective tissue-scale behaviour. The stochastic computational model simulates the interaction of auxin at a molecular scale and, by analysing the gross movement of auxin from one compartment to the next, allows us to determine auxin dynamics at the tissue scale based on the mechanistic interactions of auxin at the molecular scale. Auxin-transport speeds for the deterministic (asymptotic) and stochastic models. 95% confidence intervals are also given for the stochastic model.
The model results enable us to highlight the advantages of each approach. We solved the stochastic version of the model using a P system framework: the model is written in terms of numbers of molecules and we prescribe the probability of a molecule moving between compartments. P systems are highly intuitive, and an excellent way of engaging with a biological audience. The stochastic model generates in particular both the mean and the standard deviation of the auxin concentrations, which enables us to characterise the expected variability. We also solved the model by deriving deterministic ordinary differential equations and using asymptotic methods to obtain formulae for the auxin concentrations and transport speeds. This method requires careful analysis to determine the dominant processes on each time scale, and the resulting expressions show clearly how the model parameters affect the predicted auxin concentrations and speeds. Although numerically solving the deterministic version of the model is often the quickest method of producing a solution, stochastic P system models and asymptotic analysis can provide additional insight and information that can complement, or be an alternative, to a deterministic numerical solution. The results also highlight how the experimental set up may lead to potential discrepancies between the measured auxin velocities, and, in particular, how the measured velocities will be affected by diffusion within the agar block. Auxin speeds are generally assessed by measuring the number of auxin molecules in the collecting agar block; however, we showed that on long time scales the auxin concentration in the agar block depends on the agar-block length, and the formulae for the auxin velocity and collecting-block concentration (obtained from the asymptotic analysis) are clearly not related. We could gain further understanding of the biological implications of this result by extending the model to incorporate a more accurate representation of the stem segment, for example by modelling multiple cell files with tissue-specific active transport.
Additional file 1: Further model details. This document (PDF format) provides supporting information for the main text, and gives further details on the biological parameter estimates used in the model; the derivation of the stochastic reaction constants and the related equations that governing the deterministic model; the model checking of the stochastic computational model; and the asymptotic solution.