CMMSE: numerical analysis of a chemical targeting model

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.


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 reactiondiffusion 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 onedimensional 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 these approximations and the behaviour of the solution.

The mathematical model
Let us denote by ⊂ 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 x = (x j ) d j=1 and t ∈ [0, T ] be the spatial and temporal variables, respectively. The boundary of the domain = ∂ is assumed to be Lipschitz, and its outward unit normal vector is 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]: 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: 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: 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.
where A 0 , B 0 and C 0 are the given initial conditions and ∇ 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 ( ), E = H 1 ( ) and H = [L 2 ( )] d and denote by ·, · Y , ·, · E and ·, · H the respective scalar products in these spaces (with corresponding norms · Y , · E and · H ).
From now on, in order to simplify the writing, we do not indicate the explicit dependence of the functions on the spatial variable 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 ] → E, the concentration of the second chemical B : [0, T ] → E and the concentration of the third chemical C : In the following we state that the previous variational problem admits a unique solution.
Theorem 1 If we assume that the coefficients λ 3 , λ 4 and λ 2 , and the diffusion coefficients δ A and δ B are strictly positive, then there exists a unique solution to Problem V P with the following regularity: 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.

Numerical analysis of a fully discrete approximation
Now, we consider a fully discrete approximation of Problem VP. Firstly, we start assuming that the domain¯ is polyhedral, and denoting by T h a regular triangulation in the sense of Ciarlet [3]. Thus, we can construct the finite dimensional space E h ⊂ E as follows: where P 1 (T r) represents the space of polynomials of degree less or equal to one in element T r. 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 where 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, . . . , N . Moreover, for a continuous function f : 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 Under the conditions of Theorem 1, using the classical Lax-Milgram lemma we easily find that there exists a unique discrete solution to Problem V P 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 hk n , B n − B hk n and C n −C hk n . 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 V P hk , then it follows that there exists a positive constant M > 0, independent of the discretization parameters h and k, such that Proof First, taking as a function test v h = A hk n in discrete variational equation (11) we find that Keeping in mind that using several times Cauchy-Schwarz inequality and the following Young's inequality we easily find that Proceeding in a similar form, we obtain the following estimates for the discrete concentrations of chemicals B hk n and C hk n : Combining the previous estimates it follows that Therefore, multiplying the above estimates by k and summing up to n we obtain Finally, using a discrete version of Gronwall's inequality (see, for instance, [2]) we conclude the desired stability property. 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 ∈ E h ⊂ E and discrete variational equation (11), then we have, for all v h ∈ E h , and therefore, it follows that, for all v h ∈ E h , Now, taking into account that where we have used the well-known Cauchy-Schwarz inequality, the Young's inequality used previously, and the definition of the truncated function φ A , we conclude that 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: Combining estimates (14)-(16) it follows that Multiplying the above estimates by k and summing up to n, we have 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 V P and V P hk , then we have the following a priori error estimates, for 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).

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

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: 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 ∈ (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 calculated as 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.

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 λ 1 .

Case 1
We start with the problem corresponding with these parameters: where parameter λ 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 λ 1 .
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.
We assume again that parameter λ 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 λ 1 .
In Fig. 5 we also show the evolution of A and B with time.
In Fig. 6 we plot the steady state solutions for variables A and B, for three values of λ 1 .