# A systems-based mathematical modelling framework for investigating the effect of drugs on solid tumours

- Cong Liu
^{1}, - J Krishnan
^{1, 2}Email author and - Xiao Yun Xu
^{1}

**8**:45

https://doi.org/10.1186/1742-4682-8-45

© Liu et al; licensee BioMed Central Ltd. 2011

**Received: **14 February 2011

**Accepted: **7 December 2011

**Published: **7 December 2011

## Abstract

### Background

Elucidating the effects of drugs on solid tumours is a highly challenging multi-level problem, since this involves many complexities associated with transport and cellular response, which in turn is characterized by highly non-linear chemical signal transduction. Appropriate systems frameworks are needed to seriously address the sources of these complexities, especially from the cellular side.

### Results

We develop a skeletal modelling framework incorporating interstitial drug transport, intracellular signal processing and cell population descriptions. The descriptions aim to appropriately capture the nature of information flow. The model is deliberately formulated to start with simple intracellular descriptions so that additional features can be incorporated in a modular fashion. Two kinds of intracellular signalling modules which describe the drug effect were considered, one a monostable switch and the other a bistable switch. Analysis of our model revealed how different drug stimuli can lead to cell killing in the tumour. Interestingly both modules considered exhibited similar trends. The effects of important parameters were also studied.

### Conclusions

We have created a predictive systems platform integrating drug transport and cellular response which can be systematically augmented to include additional layers of cellular complexity. Our results indicate that intracellular signalling models which are qualitatively different can give rise to similar behaviour to simple (and typical) stimuli, and that validating intracellular descriptions must be performed with care by considering a variety of drug stimuli.

## Keywords

## Background

The need to systematically understand the complex aspects of solid tumours is evident when one considers the potentially fatal consequences which are associated with solid tumours growing unchecked. Solid tumours are a highly complex mini-universe in themselves. They are typically fed by a complex vascular network which provides blood and nutrients. This vascular network is itself more complex and irregular than vascular networks in normal tissues. The interstitium (the region of the tumour other than the vascular network) contains the tumour cells as well as the extracellular matrix. It is worth pointing out that even such a picture masks important events that occur at different time scales. For instance, a growing tumour which is not vascularized, secretes chemicals which eventually lead to its vascularization by the process of tumour-induced angiogenesis.

The complexity of the tumour environment becomes even more relevant when one attempts to evaluate systematically the effects of anti-cancer drug on tumours. Different drugs such as doxorubicin and paclitaxel have been used (and delivered in different forms) with the aim of effectively destroying tumour cells. These drugs are typically injected into the blood stream and enter the interstitium through the capillary wall. After entering the interstitium they diffuse in the interstitial space, where they may also bind to albumin or other proteins [1]. The unbound drug may be taken up by tumour cells, upon which they can act. Clearly, a number of complexities must be considered when one attempts to develop a mechanistic understanding of the effect of drug on solid tumours. These include the complex microvasculature as well as the complex structure of the interstitium [1]. Moreover, it is necessary to understand the highly non-linear nature of the cellular response in tumours and how this is affected by the tumour microenvironment [2], including both chemical and biophysical aspects.

Several attempts have been made to model mathematically the effect of drug on solid tumours [3]. These include compartmental models describing the tumour as single or discrete compartments [4, 5], transport models focusing primarily on blood flow and drug diffusion in tumours [6, 7] and pharmacokinetic and pharmacodynamic models including varying levels of description of the intracellular response. Recent computational work has begun to focus on combining interstitial transport with drug uptake by cells [8]. While all these models provide varying levels of insights, there are no models that offer a transparent systems level description of the constituent elements, with a dynamical systems basis for the description of the cellular signalling.

In this paper we take the first steps towards developing an integrative modelling framework which combines blood flow and interstitial transport, while also systematically accounting for the complexity of the relevant signal transduction in tumour cells. In the last decade, considerable amounts of interest and research activity have focussed on unravelling the intracellular and intercellular signal transduction, under the broad umbrella of Systems Biology. The approach of systems biology is especially relevant in the current context as many aspects of the cellular response to drugs in tumours are not systematically understood, and these (along with other aspects of cellular communication and interaction) may play important and unexpected roles in the actual cellular response. A better understanding of these complex interactions could also pave the way for refining treatment options. Indeed, an increasing awareness of the need for a "systems pharmacology" approach is being advocated [9], which will require a mechanistic description of tumour response, as well as intracellular biochemical interactions, with additional elements such as variability or patient specific information built into this.

Therefore, the aim of this study is to formulate a spatially distributed model to address the effects of anti-cancer drugs on tumour tissues, with a particular focus on the complexities of cellular response. This is achieved by integrating interstitial transport and cellular response, the latter presenting a considerable challenge given the complexity of cellular signalling and the high degree of nonlinearity in cellular signal processing. A unique feature of our approach is to account for the important aspects of intracellular response in a manner which is transparent and modular, so that the models can be systematically refined and augmented. This is very important as it offers the opportunity for further expansion, including incorporation of the powerful tools of systems biology in the proposed modelling framework. Our framework may thus be viewed as one which bridges the gap between pharmacodynamic type modelling and detailed mechanistic systems biology modelling, while being firmly rooted in dynamical systems approaches. Using the models presented here, analyses can then focus on understanding the effects of potentially important contributing factors, such as the nature of the intracellular signalling, the microenvironment and how this differs from tumour to tumour, the effect of intrinsic stochasticity as well as cellular variability.

## Methods

### Modelling strategy

It is recognized that the problem under consideration features a genuinely complex system, and hence any model building attempt must make choices regarding what features are central. In the work presented here, some explicit descriptions of the cellular processes are included, expanding on the modelling description of Eikenberry [8]. However, a conundrum we must face is that despite a concerted and focussed effort leading to modelling advances, many mechanistic gaps in understanding the apoptosis network (the key step in cell killing) remain. In addition, mechanistic elucidation of the role of signalling in apoptosis was performed in normal rather than tumour cells. Therefore, it is important to incorporate cellular effects at an appropriate level, so that the key features of the cellular and tissue dynamics are included without having to wait for all the relevant cellular networks to be fully understood. This informs our modelling approach: we use a combination of coarse grained descriptions of the cellular signalling dynamics, but also examine these alongside prototypical detailed models (discussed briefly) to see if there is any essential difference between the two.

The modelling framework presented here is intended to be used as a skeletal systems platform for gradual expansion to include further cellular/flow complexity. Thus data fitting is not attempted at this stage. Our model descriptions for cells are continuum based. This is to allow the additional tools available for continuum type analysis to be utilised as further complexity is added in the future (note the comments in Hinkelmann *et al*. [10] regarding discrete models).

### Model assumptions

_{C}) corresponding to the capillary radius and the outer radius (R

_{T}) for the tumour cord. To simulate drug distribution in the interstitium, it is assumed that drug is injected directly, either at the entry point of the tumour or at some other location in the body, each resulting in a specific drug concentration at the capillary wall, which is the inner boundary of the tumour cord. Different forms of drug infusion and input drug concentrations are examined.

The main assumptions of the model are (i) Cells are initially alive and present at uniform density in the tumour-cord (ii) There is minimal cell movement in the tumour cord (iii) Prior to injection of the drug, conditions are uniform inside the tumour cord (iv) Cellular variability, stochasticity, and effects of cell-to-cell interaction (mechanical and chemical) and cell-cycle are ignored.

### Extracellular drug transport

where E refers to the extracellular drug concentration (a function of radial position and time), k_{d} is the diffusion coefficient of the drug in the interstitium while c_{t} is the cell density. The two terms in the equation which involve the cell density describe the pumping out of intracellular drug concentration (denoted by I) and the cellular uptake of the extracellular drug concentration. Both these terms are described by sigmoidal saturating functions. The last two terms refer to the binding to and unbinding from proteins such as albumin, present in the interstitium. B denotes the concentration of bound extracellular drug.

where k_{db} is the diffusion coefficient of the bound drug.

which prescribes the flux F_{E} at the capillary wall as the sum of two terms. The first term is proportional to the difference between the drug concentration in the capillary S, and the extracellular concentration at the edge of the interstitium (E(r_{c})) and depends on the permeability of the capillary wall (K). The second term describes a convective flux through the capillary wall, where the transmural velocity is described by Starling's Law.

### Tumour cell density

Note that the logistic equation admits two steady states c_{t} = 0 and c_{t} = a/b, as long as a > 0. The latter steady state is the stable steady state. However, when a < 0 the only biologically relevant state is the zero steady state.

_{1}c

_{t}) and a quadratic death rate (a

_{2}c

_{t}+ bc

_{t}

^{2}). The net growth rate is given by (a

_{1}- a

_{2})c

_{t}- bc

_{t}

^{2}, which is exactly of the above form if a

_{1}> a

_{2}. This description is adopted here. A second way in which a very similar equation results is if the growth rate is a saturating function of cell density and equal to Ac

_{t}/(A

_{1}+ c

_{t}), where A is the maximal growth rate, and the death rate is given by Bc

_{t}. Here A, A

_{1}and B are all constants. Writing out the equation for c

_{t}in terms of growth and death gives

It is clear that this equation has a very similar form to Eqn (5) (except for the denominator) and in fact the qualitative property of the steady states (the dynamical systems aspect of interest) is essentially identical to that described above. Thus the logistic equation is adopted here, as it includes growth and death with saturating effects. It can be noted that if a_{1} < a_{2} Eqn (6) has only one physical steady state which is c_{t} = 0.

### Intracellular signalling

The modelling of intracellular signalling processes forms a substantial core of systems biology. Given the complexity of signalling, it is a non-trivial issue as to what the appropriate level of description should be, more so since there are many missing biological details in the signalling and many parameters unknown. At the same time there is an urgent need to start accounting for known and nontrivial aspects of the signalling, so that a predictive platform can be built for examining the effects of different cellular (intracellular or even intercellular) features. An appropriate degree of complexity is necessary for systematic evaluation and analysis of data.

As described earlier, the strategy adopted here is to start with simpler descriptions, and increase the degree of complexity gradually. In doing so, it is also important to represent the nature of the information flow correctly, and ensure that the qualitatively important features of detailed models must be accounted for. Sample detailed models in the literature are used to check the validity of the simplified models, especially for the type of drug signals encountered by cells.

In modelling the intracellular processes involving drugs, the main process of interest is drug induced apoptosis. Here, there are two key features that any model, detailed or simplified, must reflect. Firstly there must be some threshold effect present, and secondly and even more importantly, the "switch" to apoptosis must be realized in an irreversible way. Based on the systems biology literature, two types of switches are usually observed -monostable and bistable switches. Both switches are widely observed in cellular signalling generally but have very different signal transduction properties; monostable switches are completely reversible, while bistable switches can exhibit irreversibility and this has led them to be used in the modelling of irreversible processes.

Bistable processes have indeed been used in modelling irreversible cell fate decisions leading to apoptosis [13, 14]. Such modelling uses positive feedback in the caspase network (which is known to biologically exist) as a source of generating bistability. In the case of apoptosis in particular, it is not at all obvious that the irreversible fate (cell death) is necessarily or even reasonably (at the cell fate decision level) reflected as a steady state: it is very possible that an irreversible decision can lead to the triggering of critical cellular events from which there is no turning back. Thus the irreversibility could be reflected as a simple irreversible reaction, which is triggered only under very special circumstances [15]. Noting this, two models are employed here: (1) a sequential interconnection of a monostable switch and a downstream irreversible reaction effect, and (2) a bistable switch. The latter has self-contained threshold behaviour and irreversibility. It is noted that a combination of a bistable switch with an irreversible downstream reaction could also occur, with the bistable switch providing the threshold, and the downstream reaction providing the irreversible nature of the response (discussed later).

where n denotes the Hill coefficient, k_{h} an associated constant in the Hill term, and k is a parameter representing the time scale of the response. The above equation simply depicts the element R responding in a non-linear fashion to the input I. When the Hill exponent n is increased the steady state input-output response approaches a switch-like function. The above model is one concrete representation of a (monostable) switch like function: for each input value, there is only one steady state. We mention that monostable switch responses could also be realized through other means in signalling, for instance opposing enzymes acting at saturation, or the combination of stage-wise effects in a signalling cascade. However for our purposes all these models exhibit very similar input-output characteristics.

_{1}from its inactive form.

Note that in the above, under basal conditions R_{1} is very small and can only reach appreciable levels when the R level is high enough (for a substantial period of time). In our model cell death is triggered if R_{1} reaches a particular threshold R_{o}, and this is reflected at the population level in an irreversible manner. Note that for R_{1} to be above the threshold, the intracellular drug concentration must drive the upstream signal (R) above the threshold, and keep it there for a sufficient period of time.

In the above model, R represents a typical downstream intermediate element in the signalling cascade, while R_{1} represents a signal responsible for directly triggering apoptosis (and may be regarded as the output of the cascade).

It is worth also pointing out that for both the monostable and bistable modules, our signalling models are essentially minimal, but include a key intermediate step connecting input to response. This is done so that the qualitative dynamics is correctly captured, and also so that additional complexity or other factors can be appropriately incorporated.

where K_{m1} and K_{m2} are Michaelis Menten parameters. k_{fb} is a kinetic parameter which parametrizes the feedback strength. Note that the upstream signal I enzymatically acts to trigger the switch. The constants p and q serve to set the basal level and dynamic range of the module. The factor p describes the feedback present under basal conditions in the network (i.e. even in the absence of the stimulus I): even in the absence of the stimulus, the feedback results in a bistable circuit. This allows for the fact that even after the removal of a stimulus, an irreversible change in output could occur, and is exactly the way in which bistable signal processing (in the caspase network) leading to irreversible decision making in apoptosis is invoked [13].

The bistable switch has the feature that when the signal I crosses a threshold value, a sharp increase in the steady state of R is observed (see Appendix for more details). It should be noted that other variants of bistable modules examined all possess essentially similar features and input-output characteristics. In all these cases an upstream signal leads to the triggering of the bistable switch. While the equations of different modules are different in details, the key dynamical systems aspect of interest is in fact very similar (and manifests itself in very similar qualitative behaviour of these modules to signals of the kind encountered here). Here the output of this reaction drives the same activation of R_{1} just as before.

Cell death is triggered by a threshold of R_{1} being crossed. Note that the irreversibility could arise, simply by the threshold (in the bistable module) being crossed for enough time by the upstream signal. Thus even in the case where the upstream signal is transient, a different steady state of the bistable module can be attained, signifying an irreversible response. In contrast to the monostable switch, the effect of the crossing of a threshold in R_{1} on cell population density is modelled in a reversible manner. Thus in this case the irreversibility arises in the bistable switch itself: if it is permanently switched on (for instance even in a transient signal), then this will be reflected as cell death at that location in the population level.

Suitable representative parameters are employed for both switch modules. Although detailed quantitative comparison between the two models is not attempted (given their different dynamics), we seek to examine if the qualitative differences in these models are reflected at the population level. Therefore the parameter values for the two models are chosen in such a way that the threshold (of the intracellular drug concentration) for switching on is the same in both, and that they have comparable time scales. The latter is achieved by ensuring that the time to killing is the same for a reference signal value in both models.

The only remaining point is how to represent cell death at the population level. Noting that the cell population density is described in a continuous, rather than discrete form, the triggering of the intracellular threshold is incorporated by a sharply (or completely) reduced growth rate at the population level. From the dynamics of the logistic model, it is clear that if the growth rate becomes sufficiently low (i.e. a_{1} < a_{2}), eventually all the cells will perish. The same effect could also be achieved by increasing the death rate by adding an extra linear term. Both cases gave essentially similar results. Importantly, in the computation, this threshold effect is implemented in a unidirectional (irreversible) way in the monostable model but in a reversible way in the bistable model.

^{3}for extracellular, 1 ng/10

^{5}cells for intracellular), and the tumour cell density is non-dimensionalized by a typical value (10

^{6}cells/mm

^{3}). Likewise the input signal (denoted by S) is non-dimensionalized by the same representative factor (0.001 μg/mm

^{3}). Model parameters are determined as follows: based on non-dimensionalization, the size of the domain and time scales, parameters are derived either directly or indirectly from the literature, for the drug doxorubicin. The tumour growth parameters are chosen to reflect both the range of steady state, as well as the appropriate time scale. In the intracellular dynamics, the threshold to activate cell death and the time scale of the activating pathways are of special interest, and the variations of these parameters are examined here. Essential model parameters used in the present study are summarised in Tables 1 and 2. We emphasize that many of our essential conclusions are not strongly dependent on the choice of parameters.

Parameters and values used in extracellular drug transport and uptake.

Parameter | Symbol | Value | Reference |
---|---|---|---|

Free DOX diff. coefficient | k | 0.568 mm[2]/hr | [2] |

Bound DOX diff. coefficient | k | 0.032 mm[2]/hr | [2] |

Rate of transmembrane transport | V (V | 0.28 ng/(10 | [4] |

Michaelis constant for transmembrane transport | k | 0.219 μg/ml | [4] |

Michaelis constant for transmembrane transport | k | 1.37 ng/(10 | [4] |

Free DOX-albumin binding rate | k | 3000 hr- | [2] |

DOX-albumin dissociation rate | k | 1000 hr- | [2] |

Initial tumour cell density | C | 10 | [2] |

Tumour cell growth rate | a | 0.5 day | Estimated |

Saturation constant in logistic equation | b | 0.02592 mm | Estimated |

Tumour cell natural decay rate | a | 0.24 day | Estimated |

Tumour capillary radius | R | 10 μm | [2] |

Tumour cord radius | R | 120 μm | [2] |

Parameters and values used in intracellular signal transduction modules.

Parameter | Symbol | Value | Reference |
---|---|---|---|

Bistable switch | |||

Michaelis Menten constants | V | 27 hr | [6] |

Michaelis Menten constants | V | 0.459 hr | [6] |

Michaelis Menten constants | K | 100 | [6] |

Michaelis Menten constants | K | 0.01 | [6] |

Kinetic parameter mediating feedback strength | k | 2.927 hr | [6] |

Basal parameter in the bistable switch | p | 0.7 | Estimated |

Parameter mediating input regulation in the bistable switch | q | 0.3(ng/(10 | Estimated |

Monostable switch | |||

Kinetic parameter reflecting the time scale of the response | k | 0.432 hr | Estimated |

Associated constant | k | 1 ng/(10 | Estimated |

Hill coefficient | n | 10 | Estimated |

R | k | 3.6 hr | Estimated |

R | k | 0.144 hr | Estimated |

R | R | 0.9 | Estimated |

Model simulations were performed in MATLAB by suitably discretizing the diffusion in the radial direction in a standard way using finite differences [17], and solving the resulting set of ODEs using the solver ode15 s. The effect of the discretization was checked by doubling the mesh points. Further the simulation results in a few sample cases were also checked with simulations performed in COMSOL.

## Results

We present the results of the analysis of our model by considering and contrasting the results for multiple variants of the intracellular signal transduction considered. In the first case, constant levels of drug input at the capillary wall are examined, which provides important information about how the drug signal information arrives, and is processed by the tumour cells. Following this, more complex forms of drug input are analysed, specifically in the form of a single bolus, followed by multiple bolus injections. The effects of other important auxiliary parameters are also examined. A related discussion with some analytical insight is presented in the Appendix.

*et al*[18]. Figure 4(b) shows the corresponding results with the monostable switch. Here again, the same basic trend can be observed, except that there is a relatively narrow region in parameter space which leads to partial cell death. This in turn can be explained by the fact that owing to the intrinsic dynamics the bistable switch is a more sensitive differentiator between transient signals than a monostable switch. Another point worth noting is the fact that in both the monostable and bistable cases, the width of the cell death zone is a strongly non-linear function of the injection level.

## Discussion

In this paper we have taken the first steps towards a comprehensive systems-framework to investigate the effect of drugs on solid tumours. Our ultimate goal is to combine different aspects of the fluid mechanics and drug transport as well as intracellular processes to develop a comprehensive framework for elucidating the role of different factors which affect the efficacy of anti-cancer drugs. A particular impetus comes from the area of systems biology, where detailed investigations of intracellular processes have been performed in the past decade, with the promise that the fruits of such investigations will have an impact on human health and disease control.

Our approach is motivated by the desire to create an appropriate framework which is based on a dynamical systems underpinning, which can seriously tackle different aspects of the cellular complexity, and the non-linearity of cellular signal processing. We recognize that for this a predictive systems framework is needed.

Thus our immediate focus is not to fit data (see comments below).

The systems framework proposed in this study incorporates three different levels of description in a coupled manner: the dynamics of the extracellular drug, which includes details of the drug transport; the intracellular dynamics; and the population dynamics of tumour cells. Our framework incorporates basic descriptions of what we regard as the minimal elements to obtain a mechanistic description connecting drug input to cell killing, dynamically. For simplicity, a basic cylindrical-shaped tumour-cord model is analysed which contains a capillary at the centre and an interstitium in the surrounding annular space. Different modes of drug injection are simulated and model parameters corresponding to doxorubicin are used. As a first step, blood flow in the capillary vessel is not modelled explicitly (but will be included at the next stage), instead an effective drug concentration at the capillary wall is specified to represent the mode of drug injection. The extracellular drug transport in the interstitium primarily includes drug diffusion, uptake by the cell and pumping out from the cell. The latter elements depend on the local cell density. The intracellular description also incorporates the pumping in and out of drug, and minimal descriptions of the drug-triggered cell death. In this case rather than accounting for all the biochemical details, many of which are still unknown, coarse grained descriptions of signalling pathways are employed, which account for qualitatively different characteristics of the signal transduction dynamics (one involving a monostable irreversible switch, and the other a bistable switch). Our intracellular descriptions compactly encapsulate the qualitative nature of information flow in the cell.

The population level description of cells uses a basic logistic model which captures the cell birth and death in a continuum description. Since the present model only provides a skeletal description of the processes involved, the focus here is to obtain predictions regarding qualitative trends rather than quantitative values. At the same time this is very important, in order to seriously represent appropriately the complexity from the cellular side in a manner which can be systematically augmented. Overall, our modelling takes a hybrid approach which seeks to retain the predictive advantages of mechanistic models, while coarse graining intracellular descriptions to retain the relevant input-output signalling characteristics. We also avoid casting the intracellular models in concrete biochemical terms, but will do this subsequently, using this framework as a skeleton.

Two variants of intracellular description are examined to address the important question as to whether the difference between these descriptions plays an important role in the current context. Both models (monostable and bistable) are used to analyse different forms of drug input, such as sustained drug input, single and multiple bolus injections. The numerical results demonstrate that sustained drug injection results in either the entire tumour tissue being destroyed or the entire tissue surviving, depending on the concentration of drug. In the case of a single bolus, cell death occurs in a region close to the capillary wall, and the size of this region increases with an increase in the time period or concentration of the bolus. This also means that subsequent bolus injections will result in deeper penetration of the drug and hence broadening the region where cells are destroyed. Interestingly, the essential trends predicted by the monostable and bistable switches are similar, suggesting that the complexity of cellular signal processing notwithstanding, the responses to stimuli of the kind which may be encountered in realistic drug treatment are in fact basically similar. Differences could arise when one considers more complex temporal stimuli. We extended our analysis to a case where drug is released at a location away from the tumour and the drug concentration input to the tumour is determined using a pharmacokinetic model [19]. Results (not shown here) show a much smaller cell death zone, since a fair portion of the drug has been absorbed in other parts of the body. All these demonstrate that our model can be used to systematically investigate the cellular response to different forms of drug injections.

It is clear that the proposed modelling framework is the first step and thus has a number of limitations. Firstly, the most basic descriptions are used for the three levels of the modelling. Further work to include more biochemical details of the signalling network (including the detailed dynamics of the caspase network and its regulation) is needed in order to obtain a more detailed depiction of the cellular signalling. Clearly additional elements may involve either feed forward pathways, or possibly even multiple-switch type signalling, as well as redundant pathways, which may have a quantitative effect on the dynamics of the tumour. However since the two signalling models examined here, which are qualitatively completely different, produce very similar basic trends, we expect that essential conclusions arising purely from these additional variants in the signalling are likely to play a minor role. Indeed basic simulations of both the models we employ and a more detailed model of Tyson and co-workers [14] at the ODE level reveal similar characteristics of response for the types of signals considered. Therefore, the general response of more complex models is expected to indicate essentially similar trends. Based on our analysis it can be speculated that the additional layers of internal feedback/complexity are likely to be more important in affecting the tumour response to multiple bolus stimuli which are not too spaced out. It is worth pointing out that if signalling includes both a bistable switch and a downstream irreversible element (combining features of our models, with the switch effect being realized by a bistable circuit, but the irreversibility present through a downstream irreversible switch), we expect the behaviour to be very similar to the monostable switch we examined, for simpler temporal signals; in more complex temporal signals, we expect the bistable feature of signalling to play a dominant role in the signal transduction and response (analyzed subsequently). A related point to be made is that one should be very careful about claims of validating models which include cellular signalling, from scanty data.

Secondly, the present model does not include any cell variability, heterogeneity or cellular interaction. We fully recognize that inclusion of such factors can lead to partial survival of cells. Although incorporation of these factors may allow some fraction of cells to survive in a particular region, they need to be examined systematically to understand both their individual and combined effects. Noting that the complexity and specifically the non-linearity of cellular signal transduction will play an important role in conjunction with these, the entire modelling of these factors must be examined in a systematic and thorough manner, and this will be done carefully building on this framework. A similar comment can be made on cell cycle effects (which is of special relevance to analyzing the effect of drugs which target different stages of the cell cycle). Likewise the effects of factors which confer cell resistance need to be examined. All the above mentioned factors involve additional aspects of the cellular description and response. Other elements include the description of cell movement, and the effect of realistic microvascular geometries with comprehensive fluid mechanical descriptions.

Noting that multiple combinations of these effects may be present in different tumour types, it is vital to assess the roles of these factors in a systematic and predictive way, which makes clear testable predictions, and where possible to assess the role of these factors individually and together in an unambiguous fashion. This underscores our approach in dealing with the complexity of this system.

## Conclusions

Overall we have combined three levels of description to take the first steps towards a comprehensive systems description of drug effects on tumour cells. Since a particular focus in the future is to use this platform to investigate different effects of cellular complexity, a basic model structure has been adopted in the first stage. This is vital in order to build a credible platform which seriously and systematically addresses different aspects of cellular complexity. This in turn would be a necessary bridge to effectively exploit the progress in systems biology of cellular processes and bring this to bear on the problem of understanding drug delivery processes and improving their efficacy in eliminating tumour cells.

## Appendix

In this section, we briefly discuss some features of the model equations which are relevant to the work presented here. While we will not perform a detailed and exhaustive theoretical analysis, we examine some selected aspects which complement the numerical work.

_{1}which triggers the threshold. Note that in this model, once the threshold for apoptosis is triggered, the irreversibility implies that subsequent dynamics or approaching of a steady state for the monostable module variables R and R

_{1}are irrelevant. Thus it is important to also look at the dynamics. In general, for a time varying signal I(t), we have

from which R_{1}(t) can be calculated. Now the important piece of information from the dynamic variation of R_{1} is whether or not the threshold value is crossed for a given static or time-varying input I, and at what time the threshold is crossed.

As mentioned in the text, when R_{1} crosses a threshold, a significant lowering of the growth rate in the population description is implemented. This is done explicitly in an irreversible way in the monostable model simulations, and in a reversible way in the bistable model simulations. Clearly from the discussion above, we see that both models embody the fact that if the signal is above a threshold level for a sufficient period of time, this can lead to cell killing.

We now discuss some aspects of the distributed model. Looking at the logistic model for the growth rate, we see that the basal steady state is c_{t} = a/b, where a = a_{1} **-** a_{2} with a_{1} and a_{2} being the linear growth and death rates, with a_{1} **>** a_{2}. When the apoptosis threshold is crossed in the model, a_{1} is reduced to a value less than a_{2}. Thus as a result of the cellular dynamics, in both the monostable and bistable cases, a second steady state exists which is c_{t} **=** 0, which is the attained state if the threshold is crossed (note that in our monostable model if the threshold is temporally attained, then this is the only eventual state which can be reached, while in the bistable case one can say that if the threshold remains crossed at steady state, the cell density eventually becomes 0). Thus it is possible at steady state at any location for the steady state to be the basal level or to become 0. Our simulations provide ample evidence of the fact that some region gets killed completely, while others remain/regain their basal level.

which implies that the extracellular drug concentration is spatially uniform and this depends on the level of drug being pumped in. If a persistent stimulus is maintained, the drug concentration in the interstitium attains a commensurate value. If the stimulus was applied as a bolus, it means that eventually no drug is pumped in and at steady state the extracellular drug concentration is zero everywhere. However there can be regions in the interstitium where the cell density is zero, as we have seen in the simulations and from the discussion above. Thus at steady state the extracellular drug concentration will be uniform but the tumour cell density can be non-uniform.

The dynamics of the coupled system may also be studied. Without presenting this in greater detail here, we comment that the variation of the cell density is much slower than the remaining dynamics, so that the remaining dynamics can be analyzed using a perturbation analysis where the cell density variation occurs on a slower time scale.

### Supplementary results

In this subsection, we present some additional results to supplement those presented in the main text. We briefly present two sets of results, both related to bolus injections.

In the case of the single bolus injection, the effect of a single bolus on cell killing was examined, and the effect of bolus fractionation was considered as well. Here we present results related to the effect of infusion time of a bolus injection.

## Declarations

## Authors’ Affiliations

## References

- Jang SH, Wientjes MG, Lu D, Au JL: Drug delivery and transport to solid tumors. Pharm Res. 2003, 20: 1337-1350. 10.1023/A:1025785505977.View ArticlePubMedGoogle Scholar
- Trédan O, Galmarini CM, Patel K, Tannock IF: Drug resistance and the solid tumor microenvironment. J Natl Cancer Inst. 2007, 99: 1441-1454. 10.1093/jnci/djm135.View ArticlePubMedGoogle Scholar
- Liu C, Krishnan J, Stebbing J, Xu XY: Use of mathematical models to understand anticancer drug delivery and its effect on solid tumors. Pharmacogenomics. 2011, 12: 1337-1348. 10.2217/pgs.11.71.View ArticleGoogle Scholar
- El-Kareh AW, Secomb TW: A mathematical model for comparison of bolus injection, continuous infusion, and liposomal delivery of doxorubicin to tumor cells. Neoplasia. 2000, 2: 325-338. 10.1038/sj.neo.7900096.PubMed CentralView ArticlePubMedGoogle Scholar
- El-Kareh AW, Secomb TW: Two-mechanism peak concentration model for cellular pharmacodynamics of Doxorubicin. Neoplasia. 2005, 7: 705-713. 10.1593/neo.05118.PubMed CentralView ArticlePubMedGoogle Scholar
- Teo CS, Hor Keong Tan W, Lee T, Wang C-H: Transient interstitial fluid flow in brain tumors: Effect on drug delivery. Chemical Engineering Science. 2005, 60: 4803-4821. 10.1016/j.ces.2005.04.008.View ArticleGoogle Scholar
- Zhao J, Salmon H, Sarntinoranont M: Effect of heterogeneous vasculature on interstitial transport within a solid tumor. Microvasc Res. 2007, 73: 224-236. 10.1016/j.mvr.2006.12.003.View ArticlePubMedGoogle Scholar
- Eikenberry S: A tumor cord model for doxorubicin delivery and dose optimization in solid tumors. Theor Biol Med Model. 2009, 6: 16-10.1186/1742-4682-6-16.PubMed CentralView ArticlePubMedGoogle Scholar
- Yang R, Niepel M, Mitchison TK, Sorger PK: Dissecting Variability in Responses to Cancer Chemotherapy Through Systems Pharmacology. Clin Pharmacol Ther. 2010, 88: 34-38. 10.1038/clpt.2010.96.PubMed CentralView ArticlePubMedGoogle Scholar
- Hinkelmann F, Murrugarra D, Jarrah AS, Laubenbacher R: A mathematical framework for agent based models of complex biological networks. Bull Math Biol. 2011, 73: 1583-1602. 10.1007/s11538-010-9582-8.View ArticlePubMedGoogle Scholar
- Truskey GA, Yuan F, Katz DF: Transport phenomena in biological systems. 2004, Pearson Prentice HallGoogle Scholar
- Edelstein-Keshet L: Mathematical models in biology. 1988, Society for Industrial and Applied MathematicsGoogle Scholar
- Legewie S, Bluthgen N, Herzel H: Mathematical modeling identifies inhibitors of apoptosis as mediators of positive feedback and bistability. Plos Comput Biol. 2006, 2: 1061-1073.Google Scholar
- Zhang T, Brazhnik P, Tyson J: Computational analysis of dynamical responses to the intrinsic pathway of programmed cell death. Biophys J. 2009, 97: 415-434. 10.1016/j.bpj.2009.04.053.PubMed CentralView ArticlePubMedGoogle Scholar
- Albeck J, Burke J, Spencer S, Lauffenburger D, Sorger P: Modeling a snap-action, variable-delay switch controlling extrinsic cell death. PLoS Biol. 2008, 6: 2831-2852.View ArticlePubMedGoogle Scholar
- Ferrell JE, Xiong W: Bistability in cell signaling: How to make continuous processes discontinuous, and reversible processes irreversible. Chaos. 2001, 11: 227-236. 10.1063/1.1349894.View ArticlePubMedGoogle Scholar
- Thomas JW: Numerical partial differential equations. 1995, SpringerGoogle Scholar
- Seaton D, Krishnan J: A modular systems approach to elucidating the interaction of adaptive and monostable and bistable threshold modules. IET Systems Biology. 2011, 5: 81-94. 10.1049/iet-syb.2009.0061.View ArticlePubMedGoogle Scholar
- Eksborg S, Andersson M, Domellöf L, Lönroth U: A pharmacokinetic study of adriamycin and 4'epi-adriamycin after simultaneous intra-arterial liver administration. Med Oncol Tumor Pharmacother. 1986, 3: 105-110.PubMedGoogle 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.