A reaction-diffusion model of the receptor-toxin-antibody interaction

Abstract

Background

It was recently shown that the treatment effect of an antibody can be described by a consolidated parameter which includes the reaction rates of the receptor-toxin-antibody kinetics and the relative concentration of reacting species. As a result, any given value of this parameter determines an associated range of antibody kinetic properties and its relative concentration in order to achieve a desirable therapeutic effect. In the current study we generalize the existing kinetic model by explicitly taking into account the diffusion fluxes of the species.

Results

A refined model of receptor-toxin-antibody (RTA) interaction is studied numerically. The protective properties of an antibody against a given toxin are evaluated for a spherical cell placed into a toxin-antibody solution. The selection of parameters for numerical simulation approximately corresponds to the practically relevant values reported in the literature with the significant ranges in variation to allow demonstration of different regimes of intracellular transport.

Conclusions

The proposed refinement of the RTA model may become important for the consistent evaluation of protective potential of an antibody and for the estimation of the time period during which the application of this antibody becomes the most effective. It can be a useful tool for in vitro selection of potential protective antibodies for progression to in vivo evaluation.

1. Background

The successful bio-medical application of antibodies is well-documented (see [1, 2] and references therein) and there is an ever-increasing interest in the application of antibodies for a mitigation of the effect of toxins associated with various biological threats (epidemic outbreaks or malicious releases) [3–5]. With the recent progress in bio-engineering, many antibodies with different affinity parameters have been generated. For a long time the main target of antibody design has been the antibody affinity. However, according to recent results [6], affinity, on its own, is a poor predictor of protective or therapeutic potential of an antibody. In fact, the treatment effect of an antibody can be described by a consolidated parameter which includes the reaction rates of the receptor-toxin-antibody kinetics and the relative concentration of reacting species [6]. As a result, any given value of this parameter determines an associated range of antibody kinetic properties and its relative concentration in order to achieve a desirable therapeutic effect. Analytical models, similar to those reported in [6], can be a useful tool for in vitro selection of potentially protective antibodies for progression to in vivo evaluation. They can significantly reduce the cost of research and development programs by optimizing associated experimental efforts. From this perspective, extension and validation of such models becomes an important goal for biomedical modelling which is partially addressed in the current study.

There are a number of ways of refining the simple kinetic model for the Receptor-Toxin-Antibody (RTA) system proposed in [6]. The possibilities include incorporating a mechanism of receptor recycling, complex pathways for toxin internalization or multiple receptor population [7]. The focus of our study is on incorporation of the diffusion effects in the theoretical framework of RTA, i.e. enhancement of the reaction RTA model [6] with the capability to account for the diffusion fluxes of reacting species [7]. Such enhancement not only enables the application of the RTA model in more realistic setting (i.e. instead of the simplified "well-mixed" approximation [6] the reaction-diffusive RTA model can describe propagation of toxin into a single cell or into a system of cells), but also provides a high fidelity estimation of the limiting uptake rate of toxin by a cell (especially when it is limited by diffusion). More importantly, the refined model allows consistent simulation of the so-called 'window of opportunity' (period of time after exposure to toxin when the application of an antibody is the most effective). We believe the two latter parameters (the limiting uptake rate and the 'window of opportunity') can become the key parameters in the optimization study for the future antibody design.

The incorporation of diffusion fluxes into the RTA model can be implemented based on a generalization of the well-known analytical framework for ligand-receptor binding [6–10]. From a mathematical point of view, the inclusion of diffusion terms into the RTA kinetic model leads to significant complications (system of nonlinear PDEs instead of system of ODEs), which usually prevent any analytical progress and implies numerical solutions. This was the main motivation for our approach to tackle the refined RTA model. The aim of this study is to numerically evaluate the protective properties of an antibody against a given toxin in the model of a spherical cell placed into a toxin-antibody solution. We consider the problem of the RTA interaction in the most general setting, when relative concentrations of species are arbitrary and all diffusive fluxes are taken into account (toxin, antibody and associated complexes). We calculate the antibody treatment efficiency parameter under various scenarios and identify the causes of time variation of this parameter.

We also study the RTA interaction in the 'Well-Mixed Solution' (WMS) model, i.e. when the solution of a toxin, antibody, and toxin-antibody complex is assumed to be uniformly mixed and homogeneously distributed in an extracellular space. In this case all diffusion fluxes disappear and the model can be described by Ordinary Differential Equations (ODE). It is worth noting that, since in such approach receptors are still confined to the single cell surface, our model is different from the "well-mixed" model proposed in [6] where all species are homogeneously distributed over the whole space. But in the case of a low internalization rate (i.e. low toxin inflow into a cell) the governing equations of these models are of the same type.

The paper is organized as follows. In Section 3 we introduce the reaction-diffusion model for RTA. The WMS model is presented in Section 4. The results are presented in Section 5. Conclusions and summarising remarks are presented in Section 6.

2. Notation

Ω - the extracellular domain, i.e. the problem domain where species diffuse and react (i.e. toxin, antibody, and toxin-antibody complex),

S _e - the external surface of Ω,

S _c - the cell surface (inner surface of Ω),

r₀ - the concentration of receptors on the cell surface,

θ(t, x) - the the fraction of bounded receptors,

r₀θ - the concentration of the toxin-bound receptors (confined to S _c ),

r₀(1 - θ) - the concentration of free receptors,

u _T , u _A , and u _C - the concentrations of toxin, antibody, and toxin-antibody complex,

$u_T^0,u_A^0,u_C^0$- the initial concentrations,

κ _T , κ _A , and κ _C - the diffusivities of the toxin, antibody, and toxin-antibody complex,

k₁, k_-1 - the forward and reverse constants of toxin-antibody reaction rate,

k₂ and k_-2 - the forward and reverse constants of toxin and receptor binding rate,

k₃ - the rate constant of toxin internalization,

∂ _n - the outward normal derivative on S _e or S _c ,

∂ _t = ∂/∂t,

Δ - the Laplace operator,

ψ(t) - the antibody protection factor (a relative reduction of toxin inside a cell due to application of antibody).

3. Reaction-Diffusion Model for RTA Interaction

The reaction-diffusion system for the RTA interaction can be derived based on well-known results of the receptor-ligand system (law of mass action and diffusion). By including antibody into the system we arrive at the following equations

$$\left\{\begin{array}{c}{\partial }_t{u}_T=-k_1{u}_T{u}_A+k_{-1}{u}_C+{\kappa }_T\Delta u_T,x\in \Omega ,t>0,\\ u_T{|}_{S_e}=u_T^0,t>0,\\ {\partial }_n{u}_T=\frac{r_0}{{\kappa }_T}\left(-k_2\left(1-\theta \right)u_T+k_{-2}\theta \right),x\in S_c,t>0,\\ u_T{|}_{t=0}=u_T^0,x\in \Omega ,\end{array}\right.$$

(1)

$$\left\{\begin{array}{c}{\partial }_t\theta =k_2\left(1-\theta \right)u_T-k_{-2}\theta -k_3\theta ,x\in S_c,t>0,\\ \theta {|}_{t=0}=0,x\in S_c,\end{array}\right.$$

(2)

$$\left\{\begin{array}{c}{\partial }_t{u}_A=-k_1{u}_T{u}_A+k_{-1}{u}_C+{\kappa }_A\Delta u_A,x\in \Omega ,t>0,\\ u_A{|}_{S_e}=u_A^0,t>0,\\ {\partial }_n{u}_A{|}_{S_c}=0,t>0,\\ u_A{|}_{t=0}=u_A^0,x\in \Omega ,\end{array}\right.$$

(3)

$$\left\{\begin{array}{c}{\partial }_t{u}_C=k_1{u}_T{u}_A-k_{-1}{u}_C+{\kappa }_C\Delta u_C,x\in \Omega ,t>0,\\ u_C{|}_{S_e}=0,t>0,\\ {\partial }_n{u}_C{|}_{S_c}=0,t>0,\\ u_C{|}_{t=0}=0,x\in \Omega .\end{array}\right.$$

(4)

We disregard any excretion mechanism since we assume that it is nonsignificant over the time scales of interest (i.e. internalization time, time of toxin depletion etc).

The boundary conditions at the system above correspond to a case where initially the toxin and antibody are distributed homogeneously in the extracellular domain Ω. The boundary conditions on the outer boundary of the domain are assumed to be the constant concentrations of toxin and antibody and zero concentration of toxin-antibody complex. It is worth noting that in this case the gradients of u _T , u _A , u _C are nonzero at the outer surface of the domain and they provide a time-dependent influx of species into Ω (with implication no conservation law for u _T , u _A , u _C ). Indeed, in such an approach we disregard any depletion of toxin and antibody within Ω (the depletion will be taken into account in the WNS model, see below). In a practical experiment this setup can correspond to a single cell embedded into a large volume (compartment) of toxin-antibody solution, so toxin and antibody are in excess. In this context it is also worth noting that in the real biomedical scenarios the concentration of toxin is usually very low with respect to the concentration of receptor due to the high concentration of receptors on the surface of living cells and the high toxicological effect (lethal dose) of the most toxins of interest. This implies that the condition of the excess of antibody over toxin is practically relevant and are very easy to achieve (e.g. see experimental results of [11], where the concentration of ricin was about a thousand times less than the concentration of antibody), while the condition of the excess of toxin over receptor seems to be infeasible for any in vivo situation (but the latter condition still can be used in lab experiments for the model validation).

It is worth mentioning that models similar to (1)-(4) have been extensively studied in application to biouptake of pollutants by micro-organisms, cellular nutrition, heterogeneous catalysis and analytical instrumental measurements (for comprehensive review of these studies see [12–17], and references therein). Equations (1)-(4) can be presented in non-dimensional form by using scales of τ_* (time), l (length), and u_*(concentration). By substituting new variables, $x=l\bar{x}$, $t={\tau }_{*}\bar{t}$, $r_0=l{u}_{*}{\bar{r}}_0$, $u_T=u_{*}{ū}_T$, $u_A=u_{*}{ū}_A$, $u_C=u_{*}{ū}_C$, $u_{T0}=u_{*}{ū}_T^0$, $u_{A0}=u_{*}{ū}_A^0$, ${\bar{k}}_1={\tau }_{*}{u}_{*}{k}_1$, ${\bar{k}}_2={\tau }_{*}{u}_{*}{k}_2$, ${\bar{k}}_{-1}={\tau }_{*}{k}_{-1}$, ${\bar{k}}_{-2}={\tau }_{*}{k}_{-2}$, ${\bar{k}}_3={\tau }_{*}{k}_3$, ${\bar{\kappa }}_T={\tau }_{*}{\kappa }_T{l}^{-2}$, ${\bar{\kappa }}_A={\tau }_{*}{\kappa }_A{l}^{-2}$, ${\bar{\kappa }}_C={\tau }_{*}{\kappa }_C{l}^{-2}$ into (1)-(4) we can deduce the same system, but only in non-dimensional variables. Therefore, for simplicity in what follows, we treat system (1)-(4) as non-dimensional.

The main parameter of interest is the antibody protection factor (a relative reduction of toxin attached to a cell due to application of antibody). This parameter can be defined by the following expression [6]

$$\psi \left(t\right)=\frac{{\int }_{S_c}\theta {|}_{u_A^0>0}\mathsf{\text{d}}S}{{\int }_{S_c}\theta {|}_{u_A^0=0}\mathsf{\text{d}}S}.$$

(5)

By definition 0 ≤ ψ ≤ 1 with the lower values of ψ corresponding to the more profound therapeutic effect of antibody treatment.

By employing (5) it is possible to derive a simple estimation for the saturation value of parameter ψ (i.e. for the limit t → ∞). Indeed, from (1)-(4) for the steady-state limit we can write

$$\theta ={\theta }^{\mathsf{\text{sat}}}=\frac{k_2{u}_T^{\mathsf{\text{sat}}}}{k_2{u}_T^{\mathsf{\text{sat}}}+k_{-2}+k_3}=\frac{u_T^{\mathsf{\text{sat}}}}{u_T^{\mathsf{\text{sat}}}+K_2+b},$$

(6)

where $u_T^{\mathsf{\text{sat}}}$ is the saturation concentration of toxin, K₂ = k_-2/k₂, b = k₃/k₂. Then (5) leads to ψ₁ = ψ^sat where

$${\psi }^{\mathsf{\text{sat}}}=\frac{{\theta }^{\mathsf{\text{sat}}}{|}_{u_A^0>0}}{{\theta }^{\mathsf{\text{sat}}}{|}_{u_A^0=0}}.$$

(7)

So that ψ₁ can be expressed in terms of only one 'bulk' variable $u_T^{\mathsf{\text{sat}}}\ge 0$. Indeed, the value of ψ^sat can be appreciably affected by the diffusivities of species, since κ _T , κ _A , κ _C determine the saturation value $u_T^{\mathsf{\text{sat}}}$ by virtue of Eqs. (1)-(4).

4. WMS Model for RTA Interaction

The WMS model corresponds to an assumption that all species (toxin, antibody, and toxin-antibody complex) are distributed uniformly within the domain Ω. This implies no spatial gradients of concentrations, so all diffusivity terms disappear from system (1)-(4). Contrary to (1)-(4) we also assume that there are no fluxes of species across S _e , so we account for depletion of species in the cell compartment Ω (a simple yet consistent approach that accounts for the depletion effect was proposed in [17]). The process of toxin internalization (i.e. flux of toxin through the cell surface) can be modelled in this case as a given rate of toxin removal from the whole system [9]. Then the WMS model can be translated to a system of ODEs:

$$\left\{\begin{array}{c}{\stackrel{°}{u}}_T=-k_1{u}_T{u}_A+k_{-1}{u}_C-k_4{r}_0\left(k_2\left(1-\theta \right)u_T-k_{-2}\theta \right),t>0,\\ u_T{|}_{t=0}=u_T^0,\end{array}\right.$$

(8)

$$\left\{\begin{array}{c}\stackrel{°}{\theta }=k_2\left(1-\theta \right)u_T-k_{-2}\theta -k_3\theta ,t>0,\\ \theta {|}_{t=0}=0,\end{array}\right.$$

(9)

$$\left\{\begin{array}{c}{\stackrel{°}{u}}_A=-k_1{u}_T{u}_A+k_{-1}{u}_C,t>0,\\ u_A{|}_{t=0}=u_A^0,\end{array}\right.$$

(10)

$$\left\{\begin{array}{c}{\stackrel{°}{u}}_C=k_1{u}_T{u}_A-k_{-1}{u}_C,t>0,\\ u_C{|}_{t=0}=0.\end{array}\right.$$

(11)

Here a dot is placed over the variables to represent a time derivative; k₄ = S _c /V_Ω, where S _c and V_Ω are the area of cell and the extracellular volume. For instance, for a spherical cell of radius ρ _c , V_Ω is a domain between the cell and a concentric sphere of radius ρ _e > ρ _c , $V_{\Omega }=\frac{4}{3}\pi \left({\rho }_e^3-{\rho }_c^3\right)$, $S_c=4\pi {\rho }_c^2$, and $k_4=3{\rho }_c^2/\left({\rho }_e^3-{\rho }_c^3\right)$. For a simple model of cell culture (a uniformly distributed system of cells) the average density of cell distribution, n, is approximately equal to $3∕\left(4\pi {\rho }_e^3\right)$, so we can treat the 'external' scale ρ _e as the size of a compartment occupied by an individual cell in the culture. From this perspective, the dependence of ψ(ρ _e ) presented below can provide insight into the dependence of ψ on the cell packing density in the culture since ρ _e ≈ [3/(4πn)]^1/3 (see below).

The WMS model (8)-(11) is worth comparing with the model of the RTA interaction proposed in [6] (a kinetic model of uniformly distributed chemical species and cells). Despite these models being essentially different in their geometrical setting (in our case the receptors are still confined to a surface of a single cell), their governing equations become similar in the case when toxin inflow into a cell can be neglected (i.e. low internalization rate); the latter case seems to be very typical for many practical situations [7]. The WMS model (8)-(11) being a system of ODEs is much easier to analyze and solve numerically than the full RTA model (1)-(4) but indeed the WMS model cannot be used for estimating the effect of diffusivity of species on the protective properties of antibody (since it contains no diffusivity parameters).

With toxin internalization taken into account, the WMS model has only one conservation law $u_C+u_A=u_A^0$ (internalization implies that toxin is gradually taken away from the system). However, in the case of the low internalization rate we can set k₃ = 0 and also deduce an "approximate" conservation law for toxin, viz., $u_T+u_C+k_4{r}_0\theta =u_T^0$, which is similar to one used in [6]. These conservation laws significantly simplify an analytical treatment of the WMS model. For instance, from Eqs. (7) and (8)-(11) it is possible to derive an approximate analytical expression for the saturation value of protection factor ψ^sat. Actually, for the steady-state solution of system (8)-(11) without internalization rate (k₃ = 0) it is straightforward to derive the following closed equation

$$\left(1-\theta \right)\left(u_T^0-R_0\theta -\frac{\epsilon u_A^0\theta }{1+\left(\epsilon -1\right)\theta }\right)=K_2\theta ,$$

(12)

where ε = K₂/K₁, K₁ = k_-1/k₁, K₂ = k_-2/k₂, R₀ = r₀k₄ (the same equation is given in [6] for the "well mixed" model). Then the solution of this equation enables the calculation of protection factor ψ₂ = ψ^sat by means of Eq. (7).

We solve Eq. (12) numerically and compare the numerical results with the approximate analytical predictions deduced from the asymptotic solutions of Eq. (12). Some asymptotic analysis of Eq. (12) is presented in [6]. Our range of parameters corresponds to the case $R_0∕\left(\epsilon u_A^0\right)\ll 1$ and this enables derivation of the approximate formula

$${\psi }^{\mathsf{\text{sat}}}\approx {\psi }_3=\frac{F\left(u_A^0,u_T^0\right)}{F\left(0,u_T^0\right)},$$

(13)

where $F\left(x,y\right)=\left(q_1-\sqrt{q_1^2-4q_{2}y}\right)∕\left(2q_2\right)$, q₁ = K₂ + εx - (ε - 2)y and q₂ = q₁ - (εK₂ + y).

In order to verify our estimation of ψ near the saturation limit, we also solved non-steady system (8)-(11) numerically for large time and then by employing formula (7) determined function ψ₄ = ψ^sat. Table 1 shows that for the practically important cases the expressions for ψ₂, ψ₃, and ψ₄ are in the very good agreement. Table 1 also demonstrates ψ^sat for the case where internalization rate is taken into account.

Table 1 Comparison of saturation values of ψ for WMS model: ψ^sat = ψ₂ (12) and (7), ψ^sat = ψ₃ (13), and ψ^sat = ψ₄, where ψ₄ is estimated from the solution of (8)-(11) and (7) at t = 10 000 s

5. Numerical Results

We treated system (1)-(4) numerically for the spherically symmetric domain ρ ∈ [ρ _c , ρ _e ] and t > 0 with an implicit finite-difference scheme [18]. These settings constitute the standard spherical cellular model [8–10, 15]. Our selection of the values of parameters for the model (1)-(4) was motivated by the values available in the literature [11, 19–22] with the extended range to allow exploration and illustration of the various transport regimes that are possible in the RTA system. If for some parameters (i.e. diffusivity) data were not available, then we used values from similar models [7–9] and added some ranges to cater for data uncertainty and to provide sensitivity analysis. The following values were used in most calculations [7]: u_* = 6.02 · 10¹³ cm^-3, τ_* = 1 s, r₀ = 1.6 · 10⁴/S _c , where 1.6 · 10⁴ is the total number of receptors of the cell, l = 10^-2 cm, $S_c=4\pi {\rho }_c^2=4\pi \cdot 10^{-6}\mathsf{\text{ c}}{\mathsf{\text{m}}}^2$, ${\bar{r}}_0=2.115\cdot 10^{-3}$. The values of the other parameters are given in Table 2. If values of k₁, k₂, κ _A , and κ _T differ from those given in Table 2, they are specified in the legends of plots. We expect that the chosen values of parameters were representative enough to illustrate a rich variety of possible scenarios of the evolution of the RTA system and provide a reasonable estimate of timescales of the associated dynamics. The consistent match of the numerical predictions with the specific experimental results (i.e. on the ricin-neutralising antibodies [11, 19–21]) would involve some additional assumptions about the relationship between the concentration of species and observable parameters (e.g. cellular viability) and was outside of the scope of the current paper.

Table 2 Values of parameters used in calculations

The results of the numerical solutions are presented in Figures 1, 2, 3, 4, 5, 6, 7 and Tables 1, 3. As we indicated in the Background, the main purpose of our study was to estimate the effect of diffusive parameters of the species on the protective properties of an antibody. As such, most plots are presented below to illustrate this effect.

Figure 1

Effect of variation of the scale of cell compartment and toxin diffusivity on protection factor. External radius of the cell compartment ρ _e = 2 (1) and ρ _e = 5 (2), κ _T = 10^-2 (solid line), κ _T = 10^-3 (dashed line), κ _T = 10^-4 (symbols) and $u_T^0=0.5$.

Figure 2

Effect of variation of the scale of cell compartment and toxin diffusivity on protection factor. External radius of the cell compartment ρ _e = 2 (1) and ρ _e = 5 (2), κ _T = 10^-2 (solid line), κ _T = 10^-3 (dashed line), κ _T = 10^-4 (symbols) and $u_T^0=0.3$.

Figure 3

Effect of variation of the scale of cell compartment and toxin diffusivity on protection factor. External radius of the cell compartment ρ _e = 2 (1) and ρ _e = 5 (2), κ _T = 10^-2 (solid line), κ _T = 10^-3 (dashed line), κ _T = 10^-4 (symbols) and $u_T^0=0.1$.

Figure 4

Effect of the antibody diffusivity on the antibody protection factor. Antibody diffusivity κ _A = 10^-1 (1), κ _A = 10^-2 (2), κ _A = 10^-3 (3). Horizontal lines correspond to values of ψ^sat given by Eq. (7) for curves 1 and 2.

Figure 5

Effect of toxin diffusivity on antibody protection factor. Toxin diffusivity κ _T = 10^-2 (1), κ _T = 5 · 10^-3 (2), κ _T = 10^-3 (3), κ _T = 10^-4 (4). Horizontal line corresponds to value of ψ^sat given by Eq. (7) for curve 3.

Figure 6

Behavior of antibody protection function determined by WMS model for large time. Plots demonstrate convergence of ψ to saturation limit for different values of parameters k₁ and k₂ at ρ _e = 2; k₁ = 1.3 · 10^-2, k₂ : 1.25 · 10^-2 (1), 2.5 · 10^-2 (2), 5 · 10^-2 (3). Horizontal lines correspond to values of ψ^sat given by (13).

Figure 7

Behavior of antibody protection function determined by WMS model for large time. Plots demonstrate convergence of ψ to saturation limit for different values of parameters k₁ and k₂ at ρ _e = 2; k₁ = 1.3 · 10^-2, k₂ = 1.25 · 10^-2 (1), k₁ = 1.3 · 10^-2, k₂ = 2.5 · 10^-2 (2), k₁ = 1.3 · 10^-1, k₂ = 1.25 · 10^-2 (3). Horizontal lines correspond to values of ψ^sat given by (13).

Table 3 Comparison of saturation values of ψ for model (1)-(4): ψ^sat = ψ₁ and ψ^sat = ψ₅, where ψ₁ is determined by (6) and (7) with $u_T^{\mathsf{\text{sat}}}$ estimated by (1)-(4) at t = 1000 s, while ψ₅ is determined by (7) with θ^sat estimated by (1)-(4) at t = 10 000 s

To provide insight into the relation between the diffusion transport and the protective properties of an antibody in the spherical cellular model, it is convenient to employ the theoretical framework that is well-established in ecology and electrochemistry (toxin uptake by microorganisms and performance of microelectrodes) (e.g., see [15–17] and references therein). According to [15], the steady-state flux of toxin towards a spherical cell can be estimated from the following expression

$$J\left(t\right)={\Lambda }^{-1}{u}_T\left(t\right),\Lambda =\left(\frac{1}{k_3{K}_{*}}+\frac{{\rho }_c}{{\kappa }_{*}}\right),$$

(14)

where Λ is the conductance of the system (flux-concentration ratio), u _T (t) is the concentration of toxin on the outer boundary of Ω, viz. $u_T\left(t\right)=u_T^0$ for the boundary condition of constant concentration or $u_T\left(t\right)=u_T^0\mathsf{\text{exp}}\left(-t∕{\tau }_d\right)$ for the no-flux boundary condition, κ_* is the effective diffusion of the toxin, τ _d is the depletion time of toxin in the bulk, K_* = R₀/(R₀ + K₁) [6]. It can be seen that the parameter κ_* and depletion time τ _d (if the depletion of toxin is significant) become two 'aggregated' parameters that can be used to comprehensively characterize the influence of an antibody on toxin transport in the model of spherical cell.

The term κ_* /ρ _c in Eq. (14) represents the diffusive conductance and the term K_*k₃ represents the internalization conductance [15]. The ratio of the two terms is

$$L=K_{*}{k}_3{\rho }_c∕{\kappa }_{*},$$

(15)

which is called bioavalability number [15] and can be used to characterized the regime of toxin uptake by the cell [15, 16]. If L ≪ 1 the uptake flux is fully controlled by the internalization process, while in the opposite case L ≫1 it is controlled by diffusion. Note that for the case of ricin competitive binding to cell receptors and the mono-clonal antibody 2B11 the value of parameter L ≈ 10^-2, i.e. flux is mostly controlled by internalization process. Importantly, even in the case of diffusion dominated flux the transport of toxin can be characterized by a rich variety of regimes that are parameterized based on the so-called degree of lability, so these regimes correspond to the different asymptotical values of parameters κ_*, τ _d [15–17].

A detailed analysis of various regimes of diffusion controlled transport emerging in the spherical cellular model is outside the scope of the current paper, so we briefly present here only some key points that are relevant to the understanding of our numerical simulations (for details we refer the reader to [15–17]). It can be shown that the ratio p = κ_* /κ _T is always within the range 1 ≤ p ≤ ∞ with the minimal value p = 1 corresponding to the diffusion transport of toxin without presence of antibody (i.e. κ_* = κ _T ). The latter condition together with Eqs. (14) leads to a simple estimate for the long-time asymptote of the protection factor of antibody (5)

$$\psi \left(t\right)\approx {\psi }_{*}\mathsf{\text{exp}}\left(\gamma t\right),{\psi }_{*}=\frac{1+L_0}{1+L_0∕p},$$

(16)

where $\gamma =1∕{\tau }_d-1∕{\tau }_d^0,{\tau }_d^0$ is the depletion time of toxin without antibody, L₀ = K_*k₃ρ _c /κ _T .

Equation (16) is an analogue of expression (13) that accounts for diffusion effects and toxin depletion. We observe that, depending on a value of the parameter γ, the asymptotical behavior of the protection factor can be either zero (γ < 0), infinity (γ > 0) or non-zero constant (γ = 0). For the diffusion controlled flux, L₀ ≫ 1 and ψ_* = p while for the kinetically controlled regime L₀ ≪ 1 and ψ_* = 1. This implies that by changing the diffusivity of the reacting species (i.e. by introducing an antibody) it is possible (at least in theory) to control the behavior of the protection factor ψ(t). The case of the constant toxin influx (i.e. no depletion) simply corresponds to γ = 0. We observed most of these scenarios in our numerical simulations (see below).

An interesting (and not intuitively obvious) result of expression (16) is the general inequity ψ_* ≥ 1 (more precisely 1 ≤ ψ_* ≤ 1 + L₀). This means that for the diffusion-controlled scenario and for the case when toxin depletion is not significant (e.g. for γT ≪ 1), the introduction of an antibody can only increase the flux of toxin towards the cell. This result is a clear manifestation of a possible contribution of the antibody-toxin complex to the total toxin flux, described in [15–17]. Importantly, that for a case of the fast kinetics (situation when the reaction of antibody-toxin complexation is much faster than diffusion time of reacting species) the effective diffusivity κ_* is reduced to the mean diffusivity [16, 17]

$${\kappa }_{*}=\frac{u_T{\kappa }_T+u_C{\kappa }_C}{u_T+u_C},$$

(17)

that can be significantly different from the diffusivity of a toxin u _T . It is worth noting that antibody diffusivity does not appear in this expression. The "limit of mean diffusivity" for κ_* given by (17) occurs only for the system with κ _T ≠ κ _C [16, 17] (which is usually the case because of a difference in molecular weights).

Some analytical models for the calculation of the toxin depletion times ${\tau }_d^0$, τ _d have been proposed [16]. They are quite involved, and for details, we refer the reader to the original publications. The results [16] clearly demonstrate that the parameter $\gamma =1∕{\tau }_d-1∕{\tau }_d^0$ in (16) can depend on the 'external' scale ρ _e (i.e. size of the cell 'compartment') in a quite convoluted way. As was mentioned above, the scale ρ _e can be approximately related to the packing density of cells in a culture (ρ _e ≈ [3/(4πn)]^1/3), so plots ψ(ρ _e ) can be also interpreted as a simple qualitative illustration of the effect of variation in packing density n.

The plots ψ(ρ _e ) in Figures 1, 2, 3 depict the dependence of the antibody protection factor ψ on the radius of external surface ρ _e (i.e. a size of the cell compartment) and on toxin diffusivity κ _T . We calculated ψ for two values of ρ _e (ρ _e = 2 and ρ _e = 5), for two values of κ _T (κ _T = 10^-2 and κ _T = 10^-3), and three values $u_T^0$ ($u_T^0=0.5$ Figure 1, $u_T^0=0.3$ Figure 2, $u_T^0=0.1$ Figure 3).

We believe that the analytical results (16) discussed above and the numerical examples similar to those presented in Figures 1, 2, 3 may be important for either the planning of experiments (especially in cell culture) or for the correct interpretation of experimental data, since they provide a simple estimation for the amplitude of the observable effect (protection factor) and for the timescale during which this effect can occur (~ 1/γ).

The results depicted in Figures 1, 2, 3 also provide an illustrative example of the main finding of our study: the time evolution of the protection factor ψ(t) may switch from monotonic to markedly non-monotonic behavior with a variation of diffusion parameters of the RTA model. This phenomenon is in line with the theoretical framework proposed in [15–17]) and was observed frequently in our simulations.

Following the well-established application of the concept of lability to the spherical cell model [15–17], an incorporation of diffusion effects into our model enabled the simulation of a new phenomenology, which may occur in the RTA system. For instance, with system (1)-(4) we were able to model competitive behavior of the reaction and diffusion fluxes. As was mentioned above, the latter often manifests itself in rather convoluted (non-monotonic) dependence of concentration of species and their diffusion coefficients, see Figures 4 and 5. The plots in Figures 4 and 5 depict a variety of scenarios for time evolution of ψ for the different diffusivity of toxin and antibody (other parameters were the same). We can clearly see a switch from monotonic to non-monotonic behavior as we decrease diffusivity of toxin κ _T (Figure 5). The cases of non-monotonic behavior with a profound minimum of ψ(t) provide revealing examples of the practically important concept of a 'window of opportunity' discussed in the Background. Once the function ψ(t) moves far away from its minimal value, the 'blocking' effect of an antibody markedly decreases. We observe that the 'window of opportunity' is very scenario-dependent and it can be easily estimated from the plots similar to those presented in Figure 5.

The plots presented in Figures 1, 2, 3, 4, 5 illustrate detailed insights into the transport process associated with the different behavior of the protective function ψ(t) owing to the introduction of antibody into the RTA system. From these results we can see that the non-monotonic behavior is caused by the non-monotonic supply (transport) of toxin across the compartment to the cell surface; this transport can be affected by changing the diffusivity of species. The main conclusion from our numerical results is that the relative diffusivity of species can be used to control the effect of antibody treatment during a short time after the exposure to a toxin (usually a few minutes).

The plots presented in Figures 6 and 7 show the dependence of ψ on k₁ and k₂ for the WMS model and demonstrate a possible switch of ψ from monotonic to non-monotonic behavior as k₂ grows. Calculations show that ψ(t) depends very weakly on ρ _e and this an indication that the reactions in the surface layers around the cell (diffusive and reaction) provide the dominant contribution to the transport properties of the system. This result is depicted in Figures 6 and 7 where all curves are given for ρ _e = 2 and ρ _e = 5. As time increases, function ψ(t) tends to the diffusivity-dependent asymptote ψ^sat = ψ₁ for system (1)-(4) and to the reaction-dependent asymptote ψ^sat = ψ₃ for the WMS model irrespective of its short-time behavior (which is indeed controlled by (16)). The diffusive dependency of the saturation limit of the system (1)-(4) becomes evident if we recall that at the steady state, the flux of internalized toxin (i.e. flux across S _c ) should be compensated by the diffusion influx across the outer surface S _e . In order to validate this asymptotic behavior we computed a rich set of scenarios. These results are presented in Figures 4, 5, 6, 7 and in Tables 1 and 3. In Table 3 we compare ψ₁ and ψ₅ where ψ₁ is determined by (6) and (7) with $u_T^{\mathsf{\text{sat}}}$ estimated by (1)-(4) at t = 1000 s, while ψ₅ is determined by (7) with θ^sat estimated (1)-(4) at t = 10 000 s. We observe that function ψ(t) converges to an asymptotic value, but this convergence can be rather slow.

As was suggested by one of the anonymous referees, the observable strongly non-monotonic behavior of parameter ψ(t) in some of our modeling scenarios can possibly be explained by applying the concept of dynamic speciation to the formation of a toxin-antibody complex [15–17]). In the diffusion-controlled regime the dynamic speciation (i.e. the fast toxin-antibody kinetics over diffusion time) can lead to the significant contribution to the toxin flux towards the cell and (under condition κ _C < κ _T ) can even cause a 'retardation' effect [15]. After some estimations we found this hypothesis quite reasonable. For a cell size of ρ _c ≈ 10^-5 m the diffusion time is τ _κ ≈ 0.2 s for κ ≈ 1 · 10^-9 m²s^-1. The estimation for equilibration time τ _e was derived from the rigorous theoretical framework proposed in [23] for competitive binding system (application of this framework to the toxin-receptor and toxin-antibody binding can be found in [6]). Indeed the equilibration time τ _e is a strong function of the toxin concentration; it rapidly decreases as the toxin concentration increases (reacting species can faster find each other to form a complex). If as a reference point we assume that the value of parameters correspond to the scenario of binding of ricin to receptor and to the antibody then for the toxin concentration T = 10 pM the reaction time is of order of 10 s. By further increasing the toxin concentration (five times in our simulations) it appears that we approach the transition threshold from the 'inert' to the 'dynamic' complex, so the toxin-antibody complex starts contributing to the diffusion flux. A more challenging task was to identify the scenarios where this additional contribution can be appreciable, since the total flux is mainly internalization-controlled. Nevertheless, if we recall that diffusivities and concentration of species are varied by an order of hundreds, then reaching the diffusion-controlled regime in some of our simulations looks quite feasible. A more detailed interpretation of our numerical results with the concept of dynamic speciation would require additional numerical calculations (careful estimations of equilibration time of complex for each scenario) and is outside of the scope of the current study.

6. Concluding Remarks

In summary, we have refined the RTA model developed in [6] by incorporating diffusion of reacting species in the extracellular space. By solving numerically the system of nonlinear PDEs of the model we managed to simulate a rich variety of reaction-diffusion processes that may occur in the RTA system. For various combinations of parameters (rates of reactions, diffusivity and initial concentrations) we estimated the effect of antibody on the toxin penetration into a cell and expressed the effect of the antibody treatment in terms of a non-dimensional protection factor (relative reduction of toxin concentration within a cell). We demonstrated that this factor can be a significantly non-monotonic function of time and its behavior depends on an interplay between diffusive and reaction processes in the RTA model. We also examined the time evolution of the protection factor and found that it eventually tends to a diffusivity-dependent asymptotic limit, but the convergence to this asymptote may take significant time. From this perspective, the refinement of the RTA model proposed in the present study becomes important for the consistent evaluation of protective potential of an antibody and for the estimation of the time period during which the application of this antibody becomes the most effective.

The selection of the rate constants for numerical simulations was motivated by data reported in the literature [11, 19–22], with the significant ranges of variability to provide a simple sensitivity analysis for the system under consideration. The values for other constants (i.e. diffusivity) were selected based on similarity with other models [8–10]. The chosen values of parameters enable an illustrative demonstration of a rich variety of regimes of the evolution that could occur in the RTA system. These regimes are similar to ones occurring in electrochemistry and ecological studies (performance of microelectrodes and toxin uptake by microorganisms). Further validation of the proposed model with a particular set of experimental data on toxin-neutralising antibodies (e.g. [11, 21]) would require a separate study. Such a study would include an application of a data fitting algorithm that accounts for the experimental data uncertainty as well as some additional assumptions about relationships of the model predictions (concentration of species, protection factor) with the observable quantities (i.e. cellular viability). The latter assumptions may significantly affect the experimental data fit and the evaluation of predictive skills of the proposed model. We will report on such study in a separate publication.

8. Authors contributions

AS suggested a simplified model of the RTA system. VS extended the model and proposed a consistent analytical framework for the modeling and simulation of the RTA interaction. PK implemented the implicit finite-difference scheme for the numerical solution of the refined PDE system. All authors contributed to the writing of the article, and read and approved the final manuscript.