Modelling effects of internalized antibody: a simple comparative study

Background The modelling framework is proposed to study protection properties of antibodies to neutralize the effects of the plant toxin (ricin). The present study extends our previous work by including (i) the model of intracellular transport of toxin to the Endoplasmic Reticulum and (ii) the model of the internalised antibodies (when antibody is delivered directly into the cytosol). Method Simulation of the receptor-toxin-antibody interaction is implemented by solving the systems of PDEs (advection-diffusion models) or ODEs (rate models) for the underlying transport coupled with mass-action kinetics. Results As the main application of the enhanced framework we present a comparative study of two kinds (external and internalised) of antibodies. This comparison is based on calculation of the non-dimensional protection factor using the same set of parameters (geometry, binding constants, initial concentrations of species, and total initial amount of the antibody). Conclusion This research will provide a framework for consistent evaluation and comparison of different types of antibodies for toxicological applications.


Background
The plant toxin ricin made from the seeds of the castor oil plant is highly toxic to mammalian cells. It is one of the deadliest toxins known and is classified as a potential bioterror agent for which no treatment is available. Some promising results have been shown recently in the immunotherapeutic approach, i.e. application of antibodies to neutralise the effects of ricin [1][2][3][4]. With the recent progress in bio-engineering, antibodies with high affinity have been generated. The development and production of new antibodies is still an expensive process that usually includes extensive experimental studies with continuous experimental refinement of antibodies properties. Evidently, that such a retrospective evaluation of different antibodies aiming at selection of the best candidate may become very time and resource consuming.
In order to reduce this experimental burden a simple (but scientifically consistent) modelling framework has been recently proposed [5][6][7][8][9]. This framework enables extensive theoretical optimization studies to increase the protective potential of antibodies before proceeding with targeted experimental studies. http://www.tbiomed.com/content/11 /1/11 The mechanism of ricin intracellular transport involves a number of steps each with its complex phenomenology which are well-documented (see [1][2][3][4] and Refs. therein). Ricin consists of an A (RTA) and B (RTB) chain linked by a disulphide bond. RTB binds to a cell surface receptor triggering uptake and retrograde transport to the Endoplasmic Reticulum (ER). In the ER, the RTA and RTB chains are separated and the RTA is translocated across the ER membrane into the cytosol. Subsequently RTA reaches ribosome and damages the protein production machinery of the cell resulting in the cell death. In this context, the toxin concentration in the cytoplasm near ER becomes the critical quantity to estimate the toxicological impact of ricin on the cell and evaluate the protective potential of the antibody. This was a motivation to introduce a consistent quantitative characteristic for antibody comparison (see below).
For the sake of parametrisation simplicity the coarse-grained modeling framework proposed in [5][6][7][8][9] seemingly ignores these fine details of toxin binding and internalisation. In fact, it is aimed at capturing the complexity of these processes by means of a small number of 'aggregated' rate constants that can be (or have been) evaluated experimentally or numerically. Such kind of models becomes a conventional tool in pharmacological modeling (for example, see [10] and Refs. therein). From the chemical point of view the framework is similar to one well-established in electrochemistry where it is used for estimation of uptake rates of the heavy metal ions from the environment, see [11][12][13]. A practical application of the proposed models involves a numerical (or sometimes analytical) solution of a nonlinear system of PDEs (diffusion kinetics) for a given set (or range) of antibody parameters (i.e. binding rates, concentration) to infer the effect of these parameters on the protective potential of the antibody.
In the present paper we extend our previous work [5][6][7][8][9] by refining models for intracellular transport and chemical interaction of species. Motivated by experimental studies available in the literature [1][2][3][4] and possible toxicological applications we consider two scenarios of antibody delivery. In the first scenario (below refer to as Scenario I) the antibody is placed inside the cell between the ER and cell membrane. The toxin initially is delivered outside the cell (in the extracellular domain). Then it moves toward the cell and interacts with receptors on the cell membrane. Some of toxin penetrates into the cell, where it further moves toward the ER (see Figure 1a) and eventually interacts with antibody.
In the second scenario (Scenario II) the antibody is delivered outside the cell where it initially interacts with the toxin partially neutralising it. The free toxin moves toward the cell membrane where it interacts with membrane receptors. Some of toxin penetrates into the cell where it is subsequently transported toward ER (see Figure 1b).
To model biologically complex phenomenology of species transport in cell we employ the following physical mechanisms: (i) diffusion (for species in extracellular space), (ii) diffusion and advection (direction transport via microtubule network) for toxin and only diffusion for antibody and toxin-antibody complex inside the cell, (iii) mass-action kinetics (to form complexes and describe binding processes).
Although all parts of the model have been implemented with the standard numerical algorithm (an implicit finite-difference scheme [14]), an interpretation of its output and parameters tuning requires some domain knowledge (solution of nonlinear systems of PDEs). In order to reduce dimension of the parameter space of the models and ameliorate their operational application we also considered so-called simplified (or reduced) models http://www.tbiomed.com/content/11/1/11 Figure 1 The schematic diagram of receptor-toxin-antibody system: (a) -Scenario I and (b) -Scenario II. -toxin, ♦ -antibody, conglutinated and ♦ in circle -toxin-antibody complex, external sphere -cell membrane, internal sphere -ER envelope.
for each scenario. The simplified (sometimes called compartment) models are based on the 'well-mixed' assumption which states that all species have uniform concentration across the modeling domain. This assumption eliminates the a necessity to calculate the gradient driven fluxes in the models (i.e. diffusion) leading to a significant simplification (translation from PDEs to ODEs). By comparing outputs of complete and reduced models we can infer and defensively comment on the trade-off between simplicity and fidelity in modeling approach to an evaluation of each scenario.
As the main application of this framework we present a comparative study of external (conventional) and internalised antibodies. More specifically, for the same set of parameters (geometry, binding constants, initial concentration or initial amount of antibody) we calculate the non-dimensional parameter which is a ratio of reductions of toxin concentration in the cytoplasm due to introduction of two kinds of antibodies. Here u T (t, ρ n )| A∈ i , u T (t, ρ n )| A∈ e is the concentration of internalized toxin near ER for the case of internalized (i) and external antibody (e). The time evolution of this aggregated parameter δ(t) enables a consistent comparison of two scenarios of antibody treatment (i.e. either delivered extracellularly or into the cytosol). For instance, the condition δ < 1 indicates that internalised antibody performs better than conventional antibody, while at δ > 1 the conventional antibody outperforms. This comparison is the main outcome of our study. The paper is organised as follows. In Section 'The models' we introduce models for toxin-antibody interaction outside and inside the cell. Results are discussed in Section 'Numerical results'. Conclusions and summarising remarks are presented in Section 'Conclusions'.

Notation
We use the notation of papers [5][6][7]: T and A -the toxin and antibody; C = TA -the toxin-antibody complex (nontoxic); http://www.tbiomed.com/content/11/1/11 ρ -the distance to the origin; r 0 θ -the concentration of the toxin-bound receptors; r 0 (1 − θ) -the concentration of the free receptors; u A and u C -the concentrations of the antibody and toxin-antibody complex, respectively, in the domains i or e ; u T -the toxin concentration in e ∪ i ; u 0 T , u 0e A , u 0i A -the initial concentrations of the toxin and antibody in e and i , respectively; κ T , κ A , κ C -the diffusivity of the toxin, antibody, and toxin-antibody complex, respectively; k 1 , k −1 -the forward and reverse constants of the toxin-antibody reaction rate; k 2 and k −2 -the forward and reverse binding rate constants of the toxin and receptor confined on the membrane; k -the toxin internalization rate constant from e across the membrane into the cell; γ -the toxin absorption rate constant describing toxin influx into ER; v(ρ), v > 0 -the toxin advective velocity; ∂ t = ∂/∂t, = ρ −2 ∂ ∂ρ (ρ 2 ∂ ∂ρ ) -the Laplace operator.

The models
As mentioned above we study two scenarios of antibody delivery. In Scenario I the antibody is delivered inside the cell, i.e. in domain i (space between the ER and cell membrane) while in Scenario II it is delivered outside the cell (space between cell membrane and external boundary of the extracellular domain, e ). In the Scenario I toxin being initially in e (extracellular domain) moves toward the cell and interacts with receptors on the cell membrane. Subsequently some of toxin penetrates into the intracellular domain, i , where it moves toward ER and interacts with antibody. Toxin competitively reacts with antibody in the extracellular domain (Scenario II) or inside the cell (Scenario I). This 'blocking' reaction results in a reduced toxin concentration on the ER envelope and is the main bio-chemical mechanism of toxin neutralization by antibody.
In order to extend our results from a single cell model to many-cell systems we impose the no-toxin flux boundary condition on the external surface S e of the extracellular domain e (to mimic periodicity of the infinite system). This also accounts for conservation of species in the system and enables consistent simulation of depletion effects (see [13] for details).
It is well-known that binding of toxin to antibody and toxin to receptors are reversible reactions, that can be described by the equations of mass-action kinetics (see [5][6][7][8][9][10] and Refs. therein). In the context of study of processes in the Scenario I these equations http://www.tbiomed.com/content/11/1/11 have been modified to include toxin directional transport and diffusion of all species in the intracellular domain (by adding appropriate advection and diffusion terms, e.g. see [15,16]). To describe extracellular transport and toxin binding by cell receptors we model the cell as a partially absorbing sphere with distributed binding sites that follow the Langmuir adsorption model. In the intracellular domain of the Scenario I the toxin transport is modelled by the advection-diffusion equation in which an advective velocity (drift) is introduced to account for directional transport via microtubule network. For the transport of the antibody and antibody-toxin complexes we still use the pure diffusion mechanism. We also assume that antibody and toxin-antibody complexes do not internalize and never cross the cell membrane or the ER envelope (internal concentric sphere). The schematic diagram of receptor-toxin-antibody system is presented in Figure 1.
Since the process of receptor binding is very rapid (seconds) and toxin internalization is relatively slow (hours) we can effectively decoupled the model for extracellular and intracellular domains and solve them sequentially. The relative values of time scales for receptor binding and toxin internalization have also been validated retrospectively and found to be consistent with the initial assumption of domains decoupling (see below).
The Advection-Diffusion Model used to evaluate Scenario I is an extension of the models previously published [4][5][6][7][8][9]. The following new processes have been included in our simulations: (i) the intracellular traffic that captures the diffusion of all species and the main features of the toxin microtubule transport, (ii) chemical interaction of all species inside the cell (see [15] and Refs. therein), and (iii) process of the toxin transport to ER.
In order to model toxin transport towards the ER a special boundary condition was introduced on the ER envelope. The radius of that spherical envelope was about the size of the domain occupied by the ER. It is assumed that the mixed (or radiation) boundary condition for species concentrations can capture (at least phenomenologically) the complexity of toxin kinetics inside the ER [17]. Now we can formulate the equations for the Advection-Diffusion Model (refer to as the ADM1 model). In extracellular space e we have the following set of equations: Function θ determined from solution by these equations is used as the boundary condition in the equations in intracellular domain i : The set of Eqs.
(2)-(6) composes the ADM1 model. The main parameter of interest for the Scenario I is the antibody protection factor (a relative reduction of the internalized toxin due to application of antibody) defined by the expression By definition, 0 ≤ μ 1 (t) ≤ 1 with the lower values of μ 1 (t) corresponding to the more profound therapeutic effect of antibody treatment.
To simplify the ADM1 model we employ the 'well-mixed' assumption [5]. This assumption implies that all species (toxin, antibody, and toxin-antibody complex) are uniformly distributed in the calculation domain for all time, so there is no spatial gradients of concentrations. Under this condition all flux terms disappear from the equations of ADM1. The process of toxin internalization (i.e. flux of toxin through the cell surface S m and ER envelope S n ) can be modelled as an appropriate rate of toxin removal from e and i , respectively.
It is worth emphasizing that a reduction to the 'well-mixed' (or compartment) models requires fulfilment of some condition that can be formulated in terms of a smallness of the ratio = τ κ /τ R , where τ κ and τ R is a diffusion and reaction time scale in the system, respectively (diffusion can quickly restore any spatial inhomogeneity of species distribution maintaining uniform concentration). In the simplest case τ κ ∼ a 2 /κ and τ R ∼ 1/k, where a is a characteristic length scale, κ is a scale of diffusivity and k is the scale of reaction rate. This is condition for what is called diffusion-or reaction-dominated regimes in diffusion kinetics [18]. Indeed, this condition does not always hold and we do not assume that it satisfies automatically in the context of our study. Since our models are characterized by a number of spatial scales (extracellular and intracellular domains, ER envelope), hierarchy of reaction rates and diffusion coefficients, it is very difficult (or even impossible) to formulate and validate any general criteria for the feasibility of the well-mixed assumption. For this reason we take a heuristic approach in which this assumption is validated retrospectively by comparing the output of the models, viz. the full (advection-diffusion) and reduced (compartment) models.
Under the well-mixed assumption the ADM1 model reduces to two ODEs. for extracellular domain e and three ODEs for inracelluar domain: We call systems (8), (9) as the Well-Mixed Model 1 (WM1). Analogously we formulate two models used to evaluate antibody protection properties in Scenario II. For the extracellular domain e we employ the Diffusion Model proposed in [5][6][7]: with toxin concentration in i being described by the Advection-Diffusion equation We refer to systems (10), (13) as the Advection-Diffusion Model 2 (ADM2). The main parameter of interest evaluated with this model is the antibody protection factor μ 2 (t) = u T (t, ρ n ; u 0 A )| A∈ e u T (t, ρ n ; 0) . (15)  Similarly to the Scenario I we can apply the well-mixed assumption to simplify this model and reduce it to a set of ODEs in i . We call system (16) and (17) as the Well-Mixed Model 2 (WM2).
To compare the antibody protection capability in both scenarios (i.e. to infer which scenario of antibody delivery is more beneficial for a particular type of antibody) we use function δ(t) defined in Eq. (1) in which u T (t, ρ n )| A∈ i and u T (t, ρ n )| A∈ e are the toxin concentrations near ER determined for the Scenario I and II, respectively. It is evident that the following identity holds, δ(t) = μ 1 (t)/μ 2 (t).
The Advection-Diffusion Models enables the fine spatial-temporal resolution of species concentration and provide valuable insights into the phenomenology of species interaction in our system. We found them to be useful tools for understanding the evolution of aggregated parameters μ 1 , μ 2 , δ predicted by our simulations. More specifically, by means of the Advection-Diffusion Models we were able to identify a number of important regimes of toxin penetration that may have significant implications for the assessment for the antibody protection properties.
These regimes can be best described by using the 'transport phenomena' terminology which is well-established in chemical engineering, for instance see [22,23]. Initially all species are concentrated in spherical layers between the external sphere, cell membrane, and ER envelope (i.e. (ρ m , ρ e ) and (ρ n , ρ m )). This layer structure results in strong radial gradient of toxin concentration which facilitates the development of diffusion fluxes across the system (since cell membrane and ER envelope are 'penetrable'). Without antibody toxin would eventually be flushed out of the system when it reaches the ER envelope driven by the pure diffusion mechanism. The advection velocity will change toxin transport at one part of the system (from diffusion dominated to advection dominated). How much toxin can be blocked by antibody during this transition to ER depends on a number of factors. They include the initial reactants concentration (toxin and antibody), relative value of toxin advection flux, and reaction time of toxin-antibody binding. It also depends on the availability of antibody in the areas of high toxin concentration to maximize the effective toxin-antibody binding. The last condition (local availability of antibody for binding toxin) is determined by the antibody diffusivity and can be validated by inspection of the solutions of Advection-Diffusion Models (spatial-temporal outputs). This validation enables justifiable comments on the relative reduction of toxin in the output of the system (on ER envelope) and the specific values of antibody protection parameters μ 1 , μ 2 . This approach was employed in our study to identify and investigate the most 'vulnerable' regimes of toxin penetration and comment on the shapes of functions μ 1 (t), μ 2 (t). Some of these regimes are presented in examples below.
Furthermore such a 'chemical engendering' approach provides a clear path for the system optimization. Indeed, from the formal point of view our system (cell) can be cast as a chemical reactor with the only output (product), being the toxin concentration at the ER envelope. In this context, the aim of antibody treatment is to reduce this output by distributing the minimal amount of a given antibody across the system. We anticipate that in such settings the optimization problem can be tacked in reduced parameter space by using the established framework for chemical reactor design [23,24].
The results of the numerical solutions are presented in Figures 2, 3, 4, 5 and 6. The main purpose of our study was to estimate the effects of the toxin diffusivity, advective velocity, and the 'ER interface' parameter γ (rate constant of ER absorption) on the protective properties of the antibody in both scenarios. As such, most plots illustrate the effect of these parameters. In order to provide a consistent comparison of the simulation outputs we apply the same initial concentration of antibody or the same initial amount (but different concentration) of antibody. The equality of the initial amount leads to the following relation for the initial concentration http://www.tbiomed.com/content/11/1/11 which follows from the simple geometrical arguments. All calculations were performed for ρ e = 0.18, ρ m = 0.1, ρ n = 0.02 for which η = 4.871. This condition becomes especially important for the comparison of full and reduced models. Unless stated otherwise the value of all parameters (except depicted in plots) are assumed to be fixed and equal to those given above. As mentioned, our simulations reveal a number of interesting regimes of toxin propagation. In particular, in Scenario I we observed that there was a range of parameters (γ , κ T , v) for which the diffusive flux of toxin, κ T ∂ ρ u T | ρ=ρ m = (r 0 kθ − vu T )| ρ=ρ m , becomes negative near the cell membrane, i.e., it is opposite to the advective flux directed toward ER. This may reduce the total toxin flux inside the cell and provide more favorable condition for toxin-antibody binding.  In other words this phenomenon is a direct result of toxin transport being dominated by directional advection (due to velocity v). In our simulations this regime occurred for slow toxin transport characterized by parameters γ , κ T , v and the slow forward reaction rate for toxin-antibody binding. It is worth mentioning that the similar phenomenon has been reported in toxicological studies [17]. We found that the development of this 'advection  dominated' region of toxin transport with reverse diffusive flux has a positive effect on antibody protection properties since it effectively increases the efficiency toxin-antibody binding (by increasing time available for toxin-antibody reaction). We also solved problem (2)-(9) using position-dependent advective velocity v(ρ) which is zero at the membrane and then rapidly tends to a constant value as ρ → ρ n . In this case because of the boundary condition (4) 2 , the diffusive toxin flux near the cell membrane is directed to ER. results demonstrate the general trend of δ(t) increase with the increase of κ T and γ . For function μ 1 (t) (Scenario I, Eq. (7)) we found noticeable non-monotonic behavior for all κ T and γ . This implies that there is a time interval for which the protection potential of antibody in Scenario I reaches its maximum and then decays. According to Figure 2a this time interval is of order of 500 s. Figure 2b demonstrates the evolution of antibody protection factor μ 2 (t) (Scenario II, Eq. (15)) for a range of parameters κ T and γ . Although the shape of function μ 2 is very different (it is always monotonically decreasing) we can conclude that the value of μ 2 drops twice during the first time interval (about 500 s). Simulations for γ = 0.1, 0.05 and κ T = 10 −2 , 10 −3 , 10 −4 provide almost identical outputs. This leads to an important http://www.tbiomed.com/content/11/1/11 conclusion of insignificance of parameters κ T and γ for the prediction of antibody protection factors in Scenario II which is a direct consequence of equation decoupling for extracellular and intracellular domain discussed above.
The plots in Figure 3 illustrate the effect of variation of parameters κ T and γ on the comparative protective skills of the same antibody. This effect is described in terms of metric factor δ(t), Eq. (1). Our simulations generally demonstrate the decrease of parameter δ(t) as γ or κ T decrease. The simulations have been performed for the same initial amount of antibody (Eq. (18)) as well as for the same initial concentration (results not presented here). We observe that depending on value of κ T and γ there exists a time interval in which δ(t) < 1, i.e. the application of the antibody in Scenario I is more efficient than in Scenario II. We found that this interval grows as κ T or γ decreases. Overall we found that the condition δ(t) < 1 is rather sensitive to the ER absorption rate γ . This conclusion is intuitively clear since by reducing value of γ we also reduce the toxin flux through the system facilitating more effective antibody-toxin binding which, in turn, leads to more favorable output for Scenario I. Figure 4 demonstrates some universal behaviors of function δ(t) that we identified for a broad range parameters of our study. Our simulations reveal that for some intermediate time (t ∈ (2 · 10 4 ; 10 6 ) seconds) function δ(t) may approximately saturates to a constant plateau (depending on values of γ and κ), before it rapidly decays to zero. This is illustrated by two plots presented in Figure 4 (note different time scales). In particular we found that for u 0i A = u 0e A , γ = 10 −1 , κ T = 10 −2 function δ(t) varies very insignificantly (less that 10%, between 4.9 and 5.2) for the time interval (t ∈ [2·10 4 ; 5·10 5 ]) seconds, then it decreases and vanishes as t → ∞. Note, the limit δ(t) → 0 occurs at the unpractically long times, so it has more methodologically interest. We do not have a clear explanation of the specific numerical value of this saturation limit (δ(t) ≈ 5); this will require further investigation. Figure 5 shows the comparative results of the performance of the full Advection-Diffusion models (ADM1, ADM2) and the reduced models (WM1, WM2). We compare estimations of the metric parameter δ(t) provided by these models for different toxin diffusivity κ T , ER absorbtion rate γ , and initial antibody concentration. We found that for the high toxin diffusivity (κ T ≥ 10 −2 ) the outputs of full model and reduced models are almost identical. The differences between the estimations of δ(t) increases as κ T decreases and for κ T < 10 −4 the reduced models fails to provide any meaningful predictions (Figure 5b). We found that even for κ T > 10 −4 the reduced models usually poorly handle the non-monotonic behavior of function δ(t) at short times. At the long times we observed that for the broad range of parameters the full and reduced models are reasonably aligned and can yield reasonable estimation of δ as t → ∞ (within 20% accuracy), Figure 5a.
The plots in Figure 6 illustrate the effect of the toxin advective velocity, v, on the prediction of the protection factor δ(t), Eq. (1). Two different values of toxin advective velocity are used for simulations, v = 10 −3 and v = 5 · 10 −4 with γ = 0.1. The toxin diffusivity was κ T = 10 −3 (Figure 6a) and 10 −4 (Figure 6b) and the same values of the initial antibody concentration in both scenarios were used. Only two plots in these figures correspond to the well-mixed models, since they do not involve advective velocity v. In general we observed the growth of δ(t) as advective velocity increases, so Scenario II becomes more favorable. This is intuitively clear because with an increase of http://www.tbiomed.com/content/11/1/11 v the toxin neutralization time decreases and therefore Scenario I becomes less favorable. We also noticed the improved accuracy of reduced models as v increased. So, for given κ T and γ , the difference between δ(t) determined by Advection-Diffusion models and Well-Mixed ones is small only if advective velocity, v, is small. But, for small γ , it dramatically increases as v grows, while for large γ it increases but not so dramatically. This is due to the switching to the advection dominated transport of toxin in the system.

Conclusions
We present a new modeling framework to evaluate the protective potential of various antibodies. This framework is based on the equation for the diffusion-and advectiondiffusion transport of all species in the extracellular and intracellular domains and include mass-action kinetics for toxin, receptor and antibody. The advection term and radiation boundary condition on the Endoplasmic Reticulum envelope were used to model the toxin transport via the microtubule network. We estimated the protection factor of antibody (relative reduction of toxin concentration near the ER) for two scenarios of operational relevance, viz., when the antibody is delivered externally to the cell and when it is delivered directly into the cytosol. To provide a consistent comparison of antibody performance we estimated the evolution of the metric factor δ(t), which is defined as the ratio of the protections factors of an antibody for each scenarios. Based on this definition the condition δ(t) < 1 over a time span implies that the internally introduced antibody has better protective capability during this time while in the opposite case (δ(t) > 1) the internalized antibody is more effective. In this study simulated both cases. Our models reveal that depending on values of some parameters (primary toxin diffusivity and ER absorbtion rate) the plot of function δ(t) can significantly change its shape while keeping the same short-and long-time limits, by undergoing a saturation stage at the intermediate times with a universal value of δ ≈ 5.
To reduce the computational burden of the Advection-Diffusion Model (a system of nonlinear PDEs) we investigated its possible simplification. By employing the well-mixed assumption we reduced the Advection-Diffusion Model to the compartment or rate model (set of ODEs) which can be easily solved with any numerical solver. For an extended range of parameters we found that predictions of the compartment model can provide about ±20% accuracy in estimation of antibody protection potential (depending on values of some other parameters), but with decreasing toxin diffusivity limit the compartment model become a poor predictor and should not be used in this context. This finding can provide some appreciation of a possible trade-off between simplification and fidelity in toxicological modelling and can be an important criterion for selection of operational models.
We would like to emphasise that the comparative study reported in the present paper can be treated only as the first step for the application of the proposed modelling framework to the practical pharmacological studies. More extensive simulations (to provide statistically viable outcomes), data fitting and an established strategy for parameter optimization using toxicological observations will be necessary. http://www.tbiomed.com/content/11/1/11