CMMSE: numerical analysis of a chemical targeting model

Abstract

Treating specific tissues without affecting other regions is a difficult task. It is desirable to target the particular tissue where the chemical has its biological effect. To study this phenomenon computationally, in this work we numerically study a mathematical model which is written as a nonlinear system composed by three parabolic partial differential equations. The variables involved in the model are the concentration of the chemical, the concentration of the binding protein and the concentration of the chemical bound to the protein. Our aim is to propose a fully discrete approximation of this problem, using the Finite Element Method and a semi-implicit Euler scheme, in order to solve it numerically. This discrete problem is analysed, obtaining a discrete stability property and some a priori error estimates that show the algorithm converges linearly if the continuous solution is regular enough. Also, some representative examples are shown, as well as the numerical verification of the convergence.

1 Introduction

Over the last twenty years, there has been a big interest in the modelling of the target action of chemicals to particular issues where their biological action should be exerted. Some typical examples are, for instance, endogenous hormones, growth factors and prescribed drugs. In the paper [7], Zhang et al. introduced a model for the transport of IGF which suggested a way to achieve selective targeting to particular issues, in such a way that the degradation of the IGF-IGFBP complex could be able to regulate the free concentration of the IGF. Therefore, in their continuation work [5] they consider a simplified model involving only two molecules, the IGF and the IGF binding protein 3 (IGFBP3), and their small IGF-IGFBP3 binary complex. They also applied this model to the transport of a prodrug within a tumour. The basic idea of the model is that such binding proteins act as “carrier proteins”, forming IGF-IGFBP complexes, which prolong the half time of IGFs (see, e.g., [4, 6]). As it is pointed out in [5], this mechanism is biologically admissible as IGFBP-degrading proteases are capable of cleaving IGFBP into fragments that have low binding affinity for IGFs.

In the paper by Gardiner et al. the model is described by using three reaction-diffusion parabolic partial differential equations which are assumed to be coupled by several nonlinear terms. In their work, they first consider the particular case where the rate of formation of the complex within the tissue is small. In such case, for the one-dimensional setting, it is possible to calculate an analytical solution and to show that, under some conditions on the parameters, the maximum concentration of IGF is found in the centre of the tissue. In the case of the full model, a finite difference method, which is implemented in Matlab but not detailed in the paper, is applied and some numerical simulations are then presented. Therefore, in our work our aim is to continue the research started in [5, 7], by introducing a semi-explicit finite element approximation of the corresponding variational problem, providing its theoretical numerical analysis, which includes a discrete stability property and a priori error estimates, and to perform some numerical simulations which demonstrate the accuracy of this approximations and the behaviour of the solution.

2 The mathematical model

Let us denote by $\Omega \subset {\mathbb {R}}^d,\, d=1,2,3$ the spatial domain (d being the spatial dimension), and by $[0,T], \, T>0$ the time interval of interest. Let ${\varvec{x}}= (x_j)_{j=1}^d$ and $t \in [0, T]$ be the spatial and temporal variables, respectively. The boundary of the domain $\Gamma = \partial \Omega$ is assumed to be Lipschitz, and its outward unit normal vector is ${\varvec{\nu }}= (\nu _i)_{i=1}^d$.

In this section, following [5] we describe the mathematical model. In order to avoid some repetition, we only consider the dimensionless version of the constitutive equations. Therefore, let us denote by A, B and C the three chemicals arising in this reaction-diffusion problem. As it is pointed out in [5], these chemicals could refer to the concentration of IGF (A), the concentration of IGFBV3 (B) and the concentration of the IGF-IGFBP3 binary complex (C).

Let us define the source functions as follows [5]:

$$\begin{aligned}&\phi _A(A,B,C) = -\lambda _1\mu _A AB+\lambda _2\mu _A C, \\&\phi _B(A,B,C) = -\lambda _1\mu _B AB+\lambda _5 C, \\&\phi _C(A,B,C) = \lambda _1 AB. \end{aligned}$$

However, since these functions are nonlinear, in order to obtain the error estimates and also due to biological restrictions (all the concentrations must be greater than zero and bounded because the space (the tissue) is limited), we introduce the following truncation operators:

$$\begin{aligned} \begin{array}{l} R_1(A)=\left\{ \begin{array}{l}A\quad \text {if}\quad A\in [0,A^*],\\ 0\quad \text {if}\quad A<0,\\ A^* \quad \text {if}\quad A>A^*, \end{array}\right. \qquad R_2(B)=\left\{ \begin{array}{l}B\quad \text {if}\quad B\in [0,B^*],\\ 0\quad \text {if}\quad B<0,\\ B^* \quad \text {if}\quad B>B^*. \end{array}\right. \end{array} \end{aligned}$$

In the previous definitions of the truncation operators $R_i$, $i=1,2$ we have denoted by $A^*$ and $B^*$ the maximum concentrations of the chemicals A and B, respectively.

In this way, we may rewrite the above constitutive source functions as follows:

$$\begin{aligned}&\phi _A(A,B,C) = -\lambda _1\mu _A R_1(A)R_2(B)+\lambda _2\mu _A C, \\&\phi _B(A,B,C) = -\lambda _1\mu _B R_1(A)R_2(B)+\lambda _5 C, \\&\phi _C(A,B,C) = \lambda _1 R_1(A)R_2(B), \end{aligned}$$

where, making an abuse of the notation, we have used the same letters for the truncated source functions.

The resulting problem is written as follows (for any number of spatial dimensions):

Problem P. Find the concentration of the first chemical $A:{\bar{\Omega }} \times [0,T] \rightarrow {\mathbb {R}}$, the concentration of the second chemical $B:{\bar{\Omega }} \times [0,T] \rightarrow {\mathbb {R}}$ and the concentration of the third chemical $C:{\bar{\Omega }} \times [0,T] \rightarrow {\mathbb {R}}$ such that

$$\begin{aligned}&{\dot{A}} = \delta _A \, \nabla ^2 A -\lambda _3 A + \phi _A(A,B,C) \quad \text {in} \quad \Omega \times (0, T) , \end{aligned}$$

(1)

$$\begin{aligned}&{\dot{B}} = \delta _B \, \nabla ^2 B -\lambda _4B + \phi _B(A,B,C) \quad \text {in} \quad \Omega \times (0, T) , \end{aligned}$$

(2)

$$\begin{aligned}&{\dot{C}} = \nabla ^2 C -\lambda _2 C + \phi _C(A,B,C) \quad \text {in} \quad \Omega \times (0, T) , \end{aligned}$$

(3)

$$\begin{aligned}&\frac{\partial \, A}{\partial {\varvec{\nu }}}=\frac{\partial \, B}{\partial {\varvec{\nu }}}=\frac{\partial \, C}{\partial {\varvec{\nu }}} = 0 \quad \text {on} \quad \Gamma \times (0, T) , \end{aligned}$$

(4)

$$\begin{aligned}&A({\varvec{x}},0) = A_0({\varvec{x}}),\quad B({\varvec{x}},0) = B_0({\varvec{x}}) \quad \text {for a.e.} \quad {\varvec{x}}\in \Omega , \end{aligned}$$

(5)

$$\begin{aligned}&C({\varvec{x}},0) = C_0({\varvec{x}}) \quad \text {for a.e.} \quad {\varvec{x}}\in \Omega , \end{aligned}$$

(6)

where $A_0$, $B_0$ and $C_0$ are the given initial conditions and $\nabla ^2$ represents the Laplacian operator.

Now, we derive the variational formulation of the above biological problem. To this purpose, let us define the variational spaces $Y = L^2(\Omega )$, $E=H^1(\Omega )$ and ${\mathcal {H}}= [L^2(\Omega )]^d$ and denote by $\big ( \cdot , \cdot \big )_Y$, $\big ( \cdot , \cdot \big )_E$ and $\big ( \cdot , \cdot \big )_{\mathcal {H}}$ the respective scalar products in these spaces (with corresponding norms $\Vert \cdot \Vert _Y$, $\Vert \cdot \Vert _E$ and $\Vert \cdot \Vert _{\mathcal {H}}$).

From now on, in order to simplify the writing, we do not indicate the explicit dependence of the functions on the spatial variable ${\varvec{x}}$. Thus, by using Green’s formula and boundary conditions (4) the variational formulation of Problem P is written as follows.

Problem VP. Find the concentration of the first chemical $A: [0,T] \rightarrow E$, the concentration of the second chemical $B: [0,T] \rightarrow E$ and the the concentration of the third chemical $C: [0,T] \rightarrow E$ such that $A(0) = \, A_0$, $B(0) = \, B_0$, $C(0) = \, C_0$ and, for a.e. $t \in (0,T)$ and for all $v,r,z \in E$,

$$\begin{aligned}&({\dot{A}}(t) ,\, v)_Y + \delta _A (\nabla A(t) ,\, \nabla v)_{{\mathcal {H}}}+\lambda _3 (A(t) ,\, v)_{Y}\nonumber \\&\quad = (\phi _A(A(t), B(t), C(t)) ,\, v )_{Y} , \end{aligned}$$

(7)

$$\begin{aligned}&({\dot{B}}(t) ,\, r)_{Y} + \delta _B (\nabla B(t) ,\, \nabla r)_{{\mathcal {H}}}+\lambda _4 (B(t) ,\, r)_{Y}\nonumber \\&\quad = (\phi _B(A(t), B(t), C(t)) ,\, r )_{Y} , \end{aligned}$$

(8)

$$\begin{aligned}&({\dot{C}}(t) ,\, z)_{Y} + (\nabla C(t) ,\, \nabla z)_{{\mathcal {H}}}+\lambda _2 (C(t) ,\, z)_{Y}\nonumber \\&\quad = (\phi _C(A(t), B(t), C(t)) ,\, z )_{Y}. \end{aligned}$$

(9)

In the following we state that the previous variational problem admits a unique solution.

Theorem 1

If we assume that the coefficients $\lambda _3$, $\lambda _4$ and $\lambda _2$, and the diffusion coefficients $\delta _A$ and $\delta _B$ are strictly positive, then there exists a unique solution to Problem VP with the following regularity:

$$\begin{aligned} A,\,B,\,C\in C^1([0,T];Y)\cap C([0,T];E). \end{aligned}$$

The proof of the above theorem is a little bit technical and it is based on well-known results on evolution equations with monotone operators (see, e.g., [1]), and a direct application of the fixed-point theorem. However, for the sake of simplicity and since, in this paper, we focus on the numerical analysis of this problem, we skip the details of the proof.

3 Numerical analysis of a fully discrete approximation

Now, we consider a fully discrete approximation of Problem VP. Firstly, we start assuming that the domain ${\bar{\Omega }}$ is polyhedral, and denoting by ${\mathcal {T}}^h$ a regular triangulation in the sense of Ciarlet [3]. Thus, we can construct the finite dimensional space $E^h \subset E$ as follows: where $P_1(Tr)$ represents the space of polynomials of degree less or equal to one in element Tr. Therefore, the finite element space $E^h$ is composed of piecewise continuous affine functions. Here, $h> 0$ denotes the spatial discretization parameter. Moreover, we assume that the discrete initial conditions are given by

$$\begin{aligned} A_0^h = {{{\mathcal {P}}}}^h A_0,\quad B_0^h = {{{\mathcal {P}}}}^h B_0,\quad C_0^h = {{{\mathcal {P}}}}^h C_0, \end{aligned}$$

(10)

where ${{{\mathcal {P}}}}^h$ is the classical finite element interpolation operator over $E^h$ (see, e.g., [3]).

Secondly, we consider a uniform partition of the time interval [0, T] with step size $k = T / N$ and nodes $t_n = n \, k$ for $n = 0, 1, \ldots , N$. Moreover, for a continuous function $f:[0,T]\rightarrow E$, we denote $f_n=f(t_n)$.

Finally, using a hybrid combination of the implicit and the explicit Euler schemes, we obtain the fully discrete approximation of variational problem VP in the following form.

Problem VP $^{hk}$ Find the discrete concentration of the first chemical $A^{hk} = \{ A^{hk}_n\}_{n=0}^N\subset E^h$, the discrete concentration of the second chemical $B^{hk} = \{ B^{hk}_n\}_{n=0}^N\subset E^h$ and the discrete concentration of the third chemical $C^{hk} = \{ C^{hk}_n\}_{n=0}^N\subset E^h$ such that $A_0^{hk} = \, A_0^h$, $B_0^{hk} = \, B_0^h$, $C_0^{hk} = \, C_0^h$, for $n=1,\ldots , N$ and for all $v^h,r^h,z^h \in E^h$,

$$\begin{aligned}&((A^{hk}_{n}-A_{n-1}^{hk})/k ,\, v^h)_{Y} + \delta _A (\nabla A^{hk}_{n} ,\, \nabla v^h)_{{\mathcal {H}}}+\lambda _3 (A^{hk}_{n},\, v^h){Y}\nonumber \\&\quad = (\phi _A(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, v^h)_{Y} , \end{aligned}$$

(11)

$$\begin{aligned}&\big ( (B^{hk}_{n}-B_{n-1}^{hk})/k ,\, r^h \big )_{Y} + \delta _B \big ( \nabla B^{hk}_{n} ,\, \nabla r^h \big )_{{\mathcal {H}}}+\lambda _4 \big ( B^{hk}_{n},\, r^h \big )_{Y}\nonumber \\&\quad = \big ( \phi _B(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, r^h \big )_{Y} , \end{aligned}$$

(12)

$$\begin{aligned}&\big ( (C^{hk}_{n}-C_{n-1}^{hk})/k ,\, z^h \big )_{Y} + \big ( \nabla C^{hk}_{n} ,\, \nabla z^h \big )_{{\mathcal {H}}}+\lambda _2 \big ( C^{hk}_{n},\, z^h \big )_{Y}\nonumber \\&\quad = \big ( \phi _C(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, z^h \big )_{Y}. \end{aligned}$$

(13)

Under the conditions of Theorem 1, using the classical Lax-Milgram lemma we easily find that there exists a unique discrete solution to Problem $VP^{hk}$.

In the rest of the section, our aim is prove a discrete stability result and to derive some a priori error estimates on the numerical errors $A_n-A_n^{hk}$, $B_n-B_n^{hk}$ and $C_n-C_n^{hk}$.

First, we prove a discrete stability property on the discrete solution.

Theorem 2

Let the assumptions of Theorem 1 hold. If we denote by $(A^{hk},B^{hk},C^{hk})$ the solution to problem $VP^{hk}$, then it follows that there exists a positive constant $M>0$, independent of the discretization parameters h and k, such that

$$\begin{aligned} \Vert A_n^{hk}\Vert _Y^2+\Vert B_n^{hk}\Vert _Y^2+\Vert C_n^{hk}\Vert _Y^2\le M. \end{aligned}$$

Proof

First, taking as a function test $v^h=A_n^{hk}$ in discrete variational equation (11) we find that

$$\begin{aligned} \begin{array}{l} \big ( (A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n^{hk} \big )_{Y} + \delta _A \big ( \nabla A^{hk}_{n} ,\, \nabla A_n^{hk} \big )_{{\mathcal {H}}}+\lambda _3 \big ( A^{hk}_{n},\, A_n^{hk} \big )_{Y}\\ - \big ( \phi _A(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, A_n^{hk} \big )_{Y}=0. \end{array} \end{aligned}$$

Keeping in mind that

$$\begin{aligned} \begin{array}{l} \big ( (A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n^{hk} \big )_{Y} \displaystyle \ge \frac{1}{2k}\Big \{\Vert A_n^{hk}\Vert _Y^2-\Vert A_{n-1}^{hk}\Vert _Y^2\Big \},\\ \delta _A \big ( \nabla A^{hk}_{n} ,\, \nabla A_n^{hk} \big )_{{\mathcal {H}}}\ge 0,\\ \lambda _3 \big ( A^{hk}_{n},\, A_n^{hk} \big )_{Y}\ge 0, \end{array} \end{aligned}$$

using several times Cauchy-Schwarz inequality and the following Young’s inequality

$$\begin{aligned} ab\le \epsilon a^2+\frac{1}{4\epsilon }b^2, \end{aligned}$$

we easily find that

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{\Vert A_n^{hk}\Vert _Y^2-\Vert A_{n-1}^{hk}\Vert _Y^2\Big \}\le M \Big (\Vert A_n^{hk}\Vert _Y^2+ \Vert A_{n-1}^{hk}\Vert _Y^2+\Vert B_{n-1}^{hk}\Vert _Y^2+\Vert C_{n-1}^{hk}\Vert _Y^2\Big ). \end{array} \end{aligned}$$

Proceeding in a similar form, we obtain the following estimates for the discrete concentrations of chemicals $B_n^{hk}$ and $C_n^{hk}$:

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{\Vert B_n^{hk}\Vert _Y^2-\Vert B_{n-1}^{hk}\Vert _Y^2\Big \}\le M \Big (\Vert B_n^{hk}\Vert _Y^2+ \Vert A_{n-1}^{hk}\Vert _Y^2+\Vert B_{n-1}^{hk}\Vert _Y^2+\Vert C_{n-1}^{hk}\Vert _Y^2\Big ),\\ \displaystyle \frac{1}{2k}\Big \{\Vert C_n^{hk}\Vert _Y^2-\Vert C_{n-1}^{hk}\Vert _Y^2\Big \}\le M \Big (\Vert C_n^{hk}\Vert _Y^2+ \Vert A_{n-1}^{hk}\Vert _Y^2+\Vert B_{n-1}^{hk}\Vert _Y^2+\Vert C_{n-1}^{hk}\Vert _Y^2\Big ). \end{array} \end{aligned}$$

Combining the previous estimates it follows that

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{\Vert A_n^{hk}\Vert _Y^2-\Vert A_{n-1}^{hk}\Vert _Y^2\Big \}+ \frac{1}{2k}\Big \{\Vert B_n^{hk}\Vert _Y^2-\Vert B_{n-1}^{hk}\Vert _Y^2\Big \}\\ \qquad +\displaystyle \frac{1}{2k}\Big \{\Vert C_n^{hk}\Vert _Y^2-\Vert C_{n-1}^{hk}\Vert _Y^2\Big \}\\ \quad \displaystyle \le M \Big (\Vert B_n^{hk}\Vert _Y^2+ \Vert A_{n-1}^{hk}\Vert _Y^2+\Vert B_{n-1}^{hk}\Vert _Y^2+\Vert C_{n-1}^{hk}\Vert _Y^2\Big ),\\ \quad \le M \Big (\Vert C_n^{hk}\Vert _Y^2+ \Vert A_{n-1}^{hk}\Vert _Y^2+\Vert B_{n-1}^{hk}\Vert _Y^2+\Vert C_{n-1}^{hk}\Vert _Y^2\Big ). \end{array} \end{aligned}$$

Therefore, multiplying the above estimates by k and summing up to n we obtain

$$\begin{aligned} \begin{array}{l} \displaystyle \Vert A_n^{hk}\Vert _Y^2+\Vert B_n^{hk}\Vert _Y^2+\Vert C_n^{hk}\Vert _Y^2\le Mk\sum _{j=1}^n \Big (\Vert B_j^{hk}\Vert _Y^2+ \Vert A_{j}^{hk}\Vert _Y^2+\Vert C_{j}^{hk}\Vert _Y^2\Big )\\ \displaystyle \qquad \qquad + M(\Vert A_0^{h}\Vert _Y^2+\Vert B_0^{h}\Vert _Y^2+\Vert C_0^{h}\Vert _Y^2). \end{array} \end{aligned}$$

Finally, using a discrete version of Gronwall’s inequality (see, for instance, [2]) we conclude the desired stability property. $\square$

Now, we derive the error estimates for the concentration of the first chemical. Subtracting variational equation (7) at time $t=t_n$ and for all $v=v^h\in E^h\subset E$ and discrete variational equation (11), then we have, for all $v^h\in E^h$,

$$\begin{aligned} \begin{array}{l} \displaystyle \big ( {\dot{A}}_n-(A^{hk}_{n}-A_{n-1}^{hk})/k ,\, v^h \big )_{Y} + \delta _A \big ( \nabla (A_n-A^{hk}_{n}) ,\, \nabla v^h \big )_{{\mathcal {H}}}\\ \quad \displaystyle +\lambda _3 \big ( A_n^{hk},v^h \big )_{Y} - \big ( \phi _A(A_{n}, B_{n}, C_n)-\phi _A(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, v^h \big )_{Y}=0, \end{array} \end{aligned}$$

and therefore, it follows that, for all $v^h\in E^h$,

$$\begin{aligned} \begin{array}{l} \displaystyle \big ( {\dot{A}}_n-(A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n-A_n^{hk} \big )_{Y} + \delta _A \big ( \nabla (A_n-A^{hk}_{n}) ,\, \nabla (A_n-A_n^{hk}) \big )_{{\mathcal {H}}}\\ \qquad \displaystyle -\big ( \phi _A(A_{n}, B_{n}, C_n)-\phi _A(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, A_n-A_n^{hk} \big )_{Y}\\ \qquad +\lambda _3 \big ( A_n^{hk},A_n-A_n^{hk} \big )_{Y} \\ \quad = \displaystyle \big ( {\dot{A}}_n-(A^{hk}_{n}-A_{n-1}^{hk})/k ,\,A_n- v^h \big )_{Y} + \delta _A \big ( \nabla (A_n-A^{hk}_{n}) ,\, \nabla (A_n-v^h) \big )_{{\mathcal {H}}}\\ \qquad \displaystyle -\big ( \phi _A(A_{n}, B_{n}, C_n)-\phi _A(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, A_n-v^h \big )_{Y} \\ \qquad +\lambda _3 \big ( A_n^{hk},A_n-v^h \big )_{Y}. \end{array} \end{aligned}$$

Now, taking into account that

$$\begin{aligned} \begin{array}{l} \big ( {\dot{A}}_n-(A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n-A_n^{hk} \big )_{Y}\\ \quad =\big ( {\dot{A}}_n-(A_{n}-A_{n-1})/k ,\, A_n-A_n^{hk} \big )_{Y}\\ \qquad \displaystyle +\big ( (A_n-A_{n-1})/k-A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n-A_n^{hk} \big )_{Y}\\ \displaystyle \quad =\big ( {\dot{A}}_n-(A_{n}-A_{n-1})/k ,\, A_n-A_n^{hk} \big )_{Y}\\ \displaystyle \qquad +\frac{1}{2k}\Big \{\Vert A_n-A_n^{hk}\Vert _Y^2-\Vert A_{n-1}-A_{n-1}^{hk}\Vert _Y^2\Big \},\\ \delta _A (\nabla (A_n-A^{hk}_{n}) ,\, \nabla (A_n-A_n^{hk}))_{{\mathcal {H}}}\ge \delta _A \Vert \nabla (A_n-A_n^{hk})\Vert _{{\mathcal {H}}}^2,\\ \Vert (\phi _A (A_{n}, B_{n}, C_n)-\phi _A(A^{hk}_{n-1}, B^{hk}_{n-1}, C^{hk}_{n-1}) ,\, v)_{Y}\Vert \\ \quad \displaystyle \le M\Big (\Vert A_n-A_{n-1}^{hk}\Vert _Y^2+\Vert B_n-B_{n-1}^{hk}\Vert _Y^2+\Vert C_n-C_{n-1}^{hk}\Vert _Y^2+\Vert v\Vert _Y^2\Big ), \end{array} \end{aligned}$$

where we have used the well-known Cauchy-Schwarz inequality, the Cauchy inequality used previously, and the definition of the truncated function $\phi _A$, we conclude that

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{ \Vert A_n-A_n^{hk}\Vert _Y^2-\Vert A_{n-1}-A_{n-1}^{hk}\Vert _Y^2\Big \}\le M\Big (\Vert A_n-A_{n-1}^{hk}\Vert _Y^2\\ \quad +\Vert B_n-B_{n-1}^{hk}\Vert _Y^2+\Vert C_n-C_{n-1}^{hk}\Vert _Y^2+\Vert A_n-A_{n}^{hk}\Vert _Y^2+\Vert A_n-v^h\Vert _E^2\\ \quad +\big ( (A_n-A_{n-1})/k-(A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n-v^h \big )_{Y}\\ \quad + \Vert {\dot{A}}_n-(A_n-A_{n-1})/k\Vert _Y^2\Big ), \end{array} \end{aligned}$$

(14)

where, here and in what follows, $M>0$ represents a positive constant which is assumed to be independent of the discretization parameters but depending on the continuous solution.

Proceeding in a similar form, we derive the a priori error estimates for the concentration of the second and third chemicals:

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{\Vert B_n-B_n^{hk}\Vert _Y^2-\Vert B_{n-1}-B_{n-1}^{hk}\Vert _Y^2\Big \}\le M\Big (\Vert A_n-A_{n-1}^{hk}\Vert _Y^2\\ \quad +\Vert B_n-B_{n-1}^{hk}\Vert _Y^2 + \Vert {\dot{B}}_n-(B_n-B_{n-1})/k\Vert _Y^2+\Vert B_n-B_{n}^{hk}\Vert _Y^2\\ \quad +\big ( (B_n-B_{n-1})/k-(B^{hk}_{n}-B_{n-1}^{hk})/k ,\, B_n-r^h \big )_{Y}\\ \quad +\Vert C_n-C_{n-1}^{hk}\Vert _Y^2+\Vert B_n-r^h\Vert _E^2\Big ), \end{array} \end{aligned}$$

(15)

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{\Vert C_n-C_n^{hk}\Vert _Y^2-\Vert C_{n-1}-C_{n-1}^{hk}\Vert _Y^2\Big \}\le M\Big (\Vert A_n-A_{n-1}^{hk}\Vert _Y^2\\ \quad +\Vert B_n-B_{n-1}^{hk}\Vert _Y^2+ \Vert {\dot{C}}_n-(C_n-C_{n-1})/k\Vert _Y^2+\Vert C_n-C_{n}^{hk}\Vert _Y^2\\ \quad +\big ( (C_n-C_{n-1})/k-(C^{hk}_{n}-C_{n-1}^{hk})/k ,\, C_n-z^h \big )_{Y}\\ \quad +\Vert C_n-C_{n-1}^{hk}\Vert _Y^2+\Vert C_n-z^h\Vert _E^2\Big ). \end{array} \end{aligned}$$

(16)

Combining estimates (14)–(16) it follows that

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{1}{2k}\Big \{\Vert A_n-A_n^{hk}\Vert _Y^2-\Vert A_{n-1}-A_{n-1}^{hk}\Vert _Y^2\Big \}\\ \displaystyle \qquad \qquad + \frac{1}{2k}\Big \{\Vert B_n-B_n^{hk}\Vert _Y^2-\Vert B_{n-1}-B_{n-1}^{hk}\Vert _Y^2\Big \} \\ \displaystyle \qquad \qquad +\frac{1}{2k}\Big \{\Vert C_n-C_n^{hk}\Vert _Y^2- \Vert C_{n-1}-C_{n-1}^{hk}\Vert _Y^2\Big \}\\ \quad \displaystyle \le M\Big (\Vert A_n-A_{n-1}^{hk}\Vert _Y^2+\Vert B_n-B_{n-1}^{hk}\Vert _Y^2+\Vert C_n-C_{n-1}^{hk}\Vert _Y^2\\ \qquad + \Vert A_n-A_{n}^{hk}\Vert _Y^2+\Vert B_n-B_{n}^{hk}\Vert _Y^2+ \Vert {\dot{A}}_n-(A_n-A_{n-1})/k\Vert _Y^2\\ \qquad +\Vert {\dot{B}}_n-(B_n-B_{n-1})/k\Vert _Y^2+\Vert {\dot{C}}_n-(C_n-C_{n-1})/k\Vert _Y^2\\ \qquad +\Vert A_n-v^h\Vert _E^2+\Vert B_n-r^h\Vert _E^2+\Vert C_n-z^h\Vert _E^2+\Vert C_n-C_{n}^{hk}\Vert _Y^2\\ \qquad +\big ( (A_n-A_{n-1})/k-(A^{hk}_{n}-A_{n-1}^{hk})/k ,\, A_n-v^h \big )_{Y}\\ \qquad +\big ( (B_n-B_{n-1})/k-(B^{hk}_{n}-B_{n-1}^{hk})/k ,\, B_n-r^h \big )_{Y}\\ \qquad +\big ( (C_n-C_{n-1})/k-(C^{hk}_{n}-C_{n-1}^{hk})/k ,\, C_n-z^h \big )_{Y}\Big ). \end{array} \end{aligned}$$

Multiplying the above estimates by k and summing up to n, we have

$$\begin{aligned} \begin{array}{l} \displaystyle \displaystyle \Vert A_n-A_n^{hk}\Vert _Y^2+\Vert B_n-B_n^{hk}\Vert _Y^2+\Vert C_n-C_n^{hk}\Vert _Y^2\\ \quad \displaystyle \le M\sum _{j=1}^n\Big (\Vert A_j-A_{j-1}^{hk}\Vert _Y^2+\Vert B_j-B_{j-1}^{hk}\Vert _Y^2+\Vert C_j-C_{j-1}^{hk}\Vert _Y^2\\ \qquad + \Vert A_j-A_{j}^{hk}\Vert _Y^2+\Vert B_j-B_{j}^{hk}\Vert _Y^2+\Vert C_j-C_{j}^{hk}\Vert _Y^2+ \Vert {\dot{A}}_j-(A_j-A_{j-1})/k\Vert _Y^2\\ \qquad +\Vert {\dot{B}}_j-(B_j-B_{j-1})/k\Vert _Y^2+\Vert {\dot{C}}_j-(C_j-C_{j-1})/k\Vert _Y^2\\ \qquad +\Vert A_j-v_j^h\Vert _E^2+\Vert B_j-r_j^h\Vert _E^2+\Vert C_j-z_j^h\Vert _E^2\\ \qquad +\big ( (A_j-A_{j-1})/k-(A^{hk}_{j}-A_{j-1}^{hk})/k ,\, A_j-v^h_j \big )_{Y}\\ \qquad +\big ( (B_j-B_{j-1})/k-(B^{hk}_{j}-B_{j-1}^{hk})/k ,\, B_j-r^h_j \big )_{Y}\\ \qquad +\big ( (C_j-C_{j-1})/k-(C^{hk}_{j}-C_{j-1}^{hk})/k ,\, C_j-z^h_j \big )_{Y}\Big ).\\ \qquad +M\Big (\Vert A_0-A_0^{h}\Vert _Y^2+\Vert B_0-B_0^{h}\Vert _Y^2+\Vert C_0-C_0^{h}\Vert _Y^2\Big ). \end{array} \end{aligned}$$

Now, considering that

$$\begin{aligned}&k\sum _{j=1}^n\big ( (A_j-A_{j-1})/k-(A^{hk}_{j}-A_{j-1}^{hk})/k ,\, A_j-v^h_j \big )_{Y}\\&\quad \displaystyle =\sum _{j=1}^n\big ( A_j-A_{j-1}-(A^{hk}_{j}-A_{j-1}^{hk}) ,\, A_j-v^h_j \big )_{Y}\\&\quad \displaystyle =\big ( A_n-A^{hk}_{n},\,A_n-v^h_n \big )_{Y}+\big ( A_0^h-A_{0},\,A_1-v^h_1 \big )_{Y}\\&\quad +\sum _{j=1}^n\big ( A_j-A^{hk}_{j},\, A_j-v^h_j-(A_{j+1}-v^h_{j+1}) \big )_{Y}, \end{aligned}$$

and the corresponding estimates for the remaining variables B and C, applying a discrete version of Gronwall’s inequality (see [2]), we conclude the following error estimates result.

Theorem 3

Let the assumptions of Theorem 1 hold. If we denote by (A, B, C) and $(A^{hk},B^{hk},C^{hk})$ the respective solutions to problems VP and $VP^{hk}$, then we have the following a priori error estimates, for all $v^h=\{v_j^{h}\}_{j=0}^N,\, r^h=\{r_j^{h}\}_{j=0}^N,\, z^h=\{z_j^{h}\}_{j=0}^N\subset E^h$,

$$\begin{aligned}&\displaystyle \max _{0\le n\le N}\Big \{ \Vert A_n-A_n^{hk}\Vert _Y^2+\Vert B_n-B_n^{hk}\Vert _Y^2+\Vert C_n-C_n^{hk}\Vert _Y^2\Big \}\nonumber \\&\quad \displaystyle \le Mk\sum _{j=1}^N\Big ( \Vert {\dot{A}}_j-(A_j-A_{j-1})/k\Vert _Y^2\nonumber \\&\qquad +\Vert {\dot{B}}_j-(B_j-B_{j-1})/k\Vert _Y^2+\Vert {\dot{C}}_j-(C_j-C_{j-1})/k\Vert _Y^2 \nonumber \\&\qquad +\Vert C_j-z^h_j\Vert _E^2+\Vert A_j-v^h_j\Vert _E^2 +\Vert B_j-r^h_j\Vert _E^2\Big )\nonumber \\&\qquad + \frac{M}{k}\sum _{j=1}^{N-1}\Big \{\Vert A_j-v_j^h-(A_{j+1}-v_{j+1}^h)\Vert _Y^2\nonumber \\&\qquad +\Vert B_j-r_j^h-(B_{j+1}-r_{j+1}^h)\Vert _Y^2 \nonumber \\&\qquad + \Vert C_j-z_j^h-(C_{j+1}-z_{j+1}^h)\Vert _Y^2\Big \}\nonumber \\&\qquad +M\max _{0\le n\le N}\Big \{\Vert A_n-v_n^h\Vert _Y^2+\Vert B_n-r_n^h\Vert _Y^2+\Vert C_n-z_n^h\Vert _Y^2\Big \}\nonumber \\&\qquad +M\Big (\Vert A_0-A_0^{h}\Vert _Y^2+\Vert B_0-B_0^{h}\Vert _Y^2 +\Vert C_0-C_0^{h}\Vert _Y^2\Big ) , \end{aligned}$$

(17)

where $M>0$ is a positive constant assumed to be independent of the discretization parameters h and k but depending on the continuous solution.

We note that we can study the convergence order from estimates (17). As an example, we have the following result which states the linear convergence of the approximation under suitable additional regularity conditions (see [2] for details regarding the estimates of the terms which are not the usually found in the finite element approximations).

Corollary 1

If we assume that the continuous solution to Problem VP has the regularity:

$$\begin{aligned} A,\,B,\,C\in C^2([0,T];Y)\cap C([0,T];H^2(\Omega )), \end{aligned}$$

and the initial conditions have the regularity

$$\begin{aligned} A_0,\,B_0,\,C_0\in H^2(\Omega ), \end{aligned}$$

then the approximations provided by Problem $VP^{hk}$ are linearly convergent; i.e., there exists a positive constant $M>0$ such that

$$\begin{aligned} \displaystyle \max _{0\le n\le N}\Big \{ \Vert A_n-A_n^{hk}\Vert _Y+\Vert B_n-B_n^{hk}\Vert _Y+\Vert C_n-C_n^{hk}\Vert _Y\Big \} \le M(h+k). \end{aligned}$$

4 Numerical results

In this section we describe some of the numerical simulations we have performed solving a one-dimensional version of Problem P.

4.1 Numerical convergence

To check numerically the result obtained in Theorem 3 we solve Problem P for several values of discretization parameters h and k using the following parameters:

$$\begin{aligned} \begin{array}{l} \Omega =(0,1),\quad \lambda _1 = 20, \quad \lambda _2 = 5, \quad \lambda _3 = 0.5, \quad \lambda _4 = 1, \quad \lambda _5 = 1, \\ \mu _A = 1, \quad \mu _B = 1, \quad \delta _A = 1, \quad \delta _B = 1. \quad \end{array} \end{aligned}$$

The initial condition was assumed constant for all three variables and equal to one, i.e. $A(x,0) = B(x,0) = C(x,0) = 1$ for $x\in (0,1)$. The final time T was 1.

Since the problem is non-linear, an analytical solution is not easy to obtain and so, we consider as exact solution a numerical solution computed with very fine discretization parameters ($h=k=10^{-6}$). The numerical errors are therefore computed as

$$\begin{aligned} \max _{0\le n\le N}\Big \{ \Vert A_n-A_n^{hk}\Vert _Y+\Vert B_n-B_n^{hk}\Vert _Y+\Vert C_n-C_n^{hk}\Vert _Y\Big \}, \end{aligned}$$

being $A_n,\, B_n,\, C_n$ such approximated discrete solution.

The results obtained are summarized in Table 1. The numerical convergence is clearly seen. The main diagonal is plotted against $h+k$ in Fig. 1. There we can see that the convergence of the algorithm is linear, as expected.

Table 1 Numerical errors ($\times 100$)

Fig. 1

Asymptotic convergence

4.2 Stationary results

In this section we explore the one-dimensional examples presented in [5]. We remark that, in all cases, the problem evolves to a steady solution, and we study the three cases presented in the mentioned reference comparing different values of $\lambda _1$.

4.2.1 Case 1

We start with the problem corresponding with these parameters:

$$\begin{aligned} \begin{array}{l} \Omega =(0,1),\quad \lambda _2 = 5, \quad \lambda _3 = 0.5, \quad \lambda _4 = 1, \quad \lambda _5 = 1, \\ \mu _A = 1, \quad \mu _B = 1, \quad \delta _A = 1, \quad \delta _B = 1, \end{array} \end{aligned}$$

where parameter $\lambda _1$ is assumed to vary, and we use the same initial conditions as in the previous example. The discretization parameters employed are $k = 0.0001$ and $h = 0.02$. In all cases we let the solution evolve until it reaches the steady state; this happens in all cases around time $t=2.7$.

In Fig. 2 we plot the steady state solutions for variables A and C, for three values of $\lambda _1$.

Fig. 2

This figure directly correlates with Fig. 6 in [5]

Furthermore, in Fig. 3 we show the evolution of the solution at the point $x=0$ for variables A and C. We recall that this point corresponds to the center of the real domain, but since there is symmetry only half of it is simulated.

Fig. 3

Evolution of the solution at $x=0$ with time

4.2.2 Case 2

Next, we study the following case:

$$\begin{aligned} \begin{array}{l} \Omega =(0,1),\quad \lambda _2 = 1, \quad \lambda _3 = 0.5, \quad \lambda _4 = 0.5, \quad \lambda _5 = 12, \\ \mu _A = 10, \quad \mu _B = 0.2, \quad \delta _A = 0.02, \quad \delta _B = 0.1. \end{array} \end{aligned}$$

We assume again that parameter $\lambda _1$ varies and we use the same initial conditions and discretization parameters. In this case, the steady state solution is not reached until time $t=13.5$.

In Fig. 4 we plot the steady state solutions for variables A and B, for three values of $\lambda _1$.

Fig. 4

This figure directly correlates with Fig. 7 in [5]

In Fig. 5 we also show the evolution of A and B with time.

Fig. 5

Evolution of the solution at $x=0$ with time

Fig. 6

This figure directly correlates with Fig. 8 in [5]

4.2.3 Case 3

Finally, we consider the case obtained with the following data:

$$\begin{aligned} \begin{array}{l} \Omega =(0,1), \quad \lambda _2 = 1, \quad \lambda _3 = 0.5, \quad \lambda _4 = 10, \quad \lambda _5 = 12, \\ \mu _A = 10, \quad \mu _B = 0.2, \quad \delta _A = 0.02, \quad \delta _B = 0.1, \end{array} \end{aligned}$$

where, again, we assume that $\lambda _1$ varies and we use the same initial condition and discretization parameters. The steady state solution for the case with $\lambda _1=0$ is reached at time $t=13.3$; while for the case with $\lambda _1=0.01$ is reached at $t=16.4$.

In Fig. 6 we plot the steady state solutions for variables A and B, for three values of $\lambda _1$.