Optimal control of drug delivery to brain tumors for a distributed parameters model

Siddhartha P. Chakrabarty Department of Mathematics, Statistics and Computer Science University of Illinois at Chicago 851 S Morgan St.; m/c 249 Chicago, IL 60607–7045, USA E–Mail: pratim@math.uic.edu URL: http://www.math.uic.edu/∼pratim/  and  Floyd B. Hanson Laboratory for Advanced Computing University of Illinois at Chicago 851 S Morgan St.; m/c 249 Chicago, IL 60607–7045, USA E–Mail: hanson@math.uic.edu URL: http://www.math.uic.edu/∼hanson/

Abstract— The growth and treatment of brain tumors is mathematically examined using a distributed parameters model. The model is a system of three coupled reaction diffusion equations involving the tumor cells, normal tissue and the drug concentration. An optimal control problem is designed, with the drug delivery rate as the control and solved to obtain the state and co-state equations as well as the regular control using a modified double shot forwardbackward method. This gives rise to a coupled system of equations with a forward state equation and a backward costate equation, which is solved using a double shot method. A numerical procedure based upon the Crank-Nicolson method is used to solve the coupled system of two three-dimension partial differential equations.

I. Introduction

The growth and control of brain tumors have been the subject of medical and scientific scrutiny for a very long time. Simply speaking a tumor, like most cancerous cells originate from a cell, that proliferates and effects its neighboring normal tissues. As the tumor cells become malignant they become more dangerous for the host. Understanding the mechanism of tumor progression is necessary for its diagnosis and treatment. As mentioned earlier, brain tumors have been studied for centuries, dating back to as early as 2500 BC [11]. The most common and deadly form of brain tumor are the gliomas, which account for more than half of the brain tumor cases. Gliomas are highly invasive and severely infiltrate the surrounding tissues [12]. Despite improved diagnostic procedures such as computerized tomography (CT) scan and magnetic resonance imaging (MRI), the benefits of such modern accessories have been restricted by the treatment options available. One major problem of administering the drugs to the brain tumor site is the blood brain barrier (BBB) [1], which exists in the human brain as a protection for the brain cells and as a restriction on the transport of water soluble substances between the blood and the central nervous system. Another problem that arises is the resection of a tumor, after the core

Work supported in part by the National Science Foundation Grants DMS–99–73231 and DMS–02–07081. The first author was also supported by the Provost’s award for Graduate Research awarded by the Graduate College at the University of Illinois at Chicago. Submitted to Proceedings of American Control Conference 2005.

mass of the tumor has been surgically removed. One way to deliver the drug is to use the drug conjugated with a nontoxic polymer like the BCNU and place it in the brain cavity for controlled release. A more traditional way would be to deliver the drug by an optimal distribution of the drug about the original tumor site. Wang et al. [13], [14] have studied the drug delivery method to tumors in three dimension for drugs like IgG and BCNU.

The paper will mostly focus on the control for the optimal distribution of the drug about the original tumor site. Before that we need to look at the mechanisms behind the growth of tumor cells as well as normal tissues and the drug concentration in the tumor. Unlike a lot of other tumors, gliomas can be highly diffusive [11]. Gatenby et al. [2] and Mansuri [9] study the mechanism of reaction diffusion in the growth of tumors. They also take into account the effects of competition for resources between the cancerous cells and the healthy tissues. Westman et al.

[15] look at the various types of tumor growth, namely exponential, logistic and Gompertz. Murray’s books [10],

[11] are excellent references for the study of different types of growth mechanisms. Also, Woodward et al. [16] study a model of glioma growth and the effects of surgical resection. In this paper we set up a fairly generalized distributed parameters model for the PDE driven system, define an objective functional to minimize drug delivery and tumor burden costs, and use a modified Lagrange multiplier method [3] for including constraints in the control problem. Finally a double shot, forward and backward, numerical method is given to solve for the state, co-state and regular control. To solve for the PDEs we use a Crank-Nicolson predictor-corrector method developed by Hanson et al. [5] and Hanson [6], for stochastic dynamic programming.

II. Mathematical Model

Let Y_{1} = n_{1}(x, t) be the density of tumor cells, Y_{2} = n_{2}(x, t) be the density of normal tissue and Y_{3} = c(x, t) be the drug concentration at any vector position x and time t. This spatio-temporal model is a system of three coupled reaction-diffusion equations.

A. Tumor Cells

It is assumed that the density of tumor cells, n_{1} = n_{1}(x, t), satisfy a reaction-diffusion equation subject to competition with the normal cells, n_{2} = n_{2}(x, t),

\frac{∂n_{1}}{∂t}=D_{1}∇_{x}^{2}[n_{1}]+a_{1}f_{1}(n_{1})−(α_{1,2}n_{2}+b_{1}(c)) n_{1}, (1)

where the tumor diffusivity is D_{1}. Let the term a_{1}f_{1}(n_{1}) be the growth rate of the tumor cells, where f_{1} could be exponential, logistic or Gompertz growth ,

f_{1}(n_{1}) = n_{1}, n_{1} (1 − \frac{n_{1}}{k_{1}} ) or n_{1} ln (\frac{k_{1}}{n_{1}} ) ,

respectively, where a_{1} is the tumor cell intrinsic growth rate and k_{1} is the tumor cell carrying capacity. Let α_{1,2} denote the death rate of the tumor cells due to competition for resources with the normal tissue. Let b_{1}(c) be the death rate of tumor cells due to drug treatment, which could be a function of the localized drug concentration c(x, t) at the tumor site, e.g., purely linear function,

b_{1}(c(x, t))=κ_{1,3}·c(x, t)

respectively.

B. Normal Tissue

Similar assumptions are made for the density of normal cells n_{2} = n_{2}(x, t) with similar coefficients. Thus, the reaction-diffusion equation for normal tissue evolution is as follows,

\frac{∂n_{2}}{∂t}=D_{2}∇_{x}^{2}[n_{2}]+a_{2}f_{2}(n_{2})−(α_{2,1}n_{1}+b_{2}(c)) n_{2}, (2)

where a_{2} is the normal cell intrinsic growth rate and the normal cell growth function f_{2} is either exponential, logistic or Gompertz growth ,

f_{2}(n_{2}) = n_{2}, n_{2} (1 − \frac{n_{2}}{k_{2}} ) or n_{2} ln (\frac{k_{2}}{n_{2}} ) ,

respectively, where k_{2} is the normal tissue carrying capacity. Note that the b_{2}(c) term indicates that some normal tissues could die as a result of the treatment and b_{2}(c) must also depend on the local drug concentration c(x, t), e.g.,

b_{2}(c(x, t)) = κ_{2,3} · c(x, t).

C. Concentration

It is assumed that the drug shows a diffusive behavior and that there is a reabsorption at the rate a_{3}. Also let u = u(x, t) be the rate at which the drug is being injected. The choice of the letter u indicates that we will use it as the control when dealing with the optimal control system. The equation for drug concentration at position x and time t is,

\frac{∂c}{∂t}=D_{3}∇_{x}^{2}[c]+a_{3}f_{3}(c)+u,  (3)

where f_{3}(c) = −c is the reabsorption function[14].

D. Global State Vector Let the global state vector be

Y(x, t) = [Y_{i}(x, t)]_{3×1} = [\begin{matrix}n_{1}(x, t) \\ n_{2}(x, t) \\ c(x, t) \end{matrix}] ,  (4)

at position x in the interior Ω of the state domain and time t on [0, t_{f}].

E. Initial and Boundary Conditions Let the initial conditions for the state be

Y(x, 0) = [\begin{matrix}n_{1}(x, 0) \\ n_{2}(x, 0) \\ c(x, 0) \end{matrix}] = [\begin{matrix}n_{1}_{0}(x) \\ n_{2}_{0}(x) \\ c_{0}(x) \end{matrix}] ≡ Y_{0}(x),  (5)

for x in Ω. Murray [11] recommends using Gaussian distribution for the initial distributions of tumors. The no flux boundary conditions are

−D(N̂ · ∇_{x})[Y](x, t)=[\begin{matrix}−D_{1}(N̂·∇_{x})[n_{1}] \\ −D_{2}(N̂·∇_{x})[n_{2}] \\ −D_{3}(N̂·∇_{x})[c] \end{matrix}](x, t)=0, (6)

for x ∈ Γ = ∂Ω, i.e., on the boundary of the domain, and for t ∈ [0, t_{f}], assuming D_{i} ≠ 0 or else the D_{i} would not be used in the condition, where N̂(x, t) is the normal to the boundary, −D_{i}∇_{x}[Y_{i}] is the flux of the ith component and the diffusion matrix,

D(x) = [\begin{matrix}D_{1} & 0 & 0 \\ 0 & D_{2} & 0 \\ 0 & 0 & D_{3} \end{matrix}] = [D_{i}δ_{i,j}]_{3×1}

is diagonal and inhomogeneous depending on the brain matter [12], where δ_{i,j} is the Kronecker delta. Note that the no flux condition at the boundary is motivated by the physical reality that the brain is a finite and closed domain.

III. Optimal Control Problem A possible objective functional is the quadratic form of

running and terminal costs,

J(u) = \frac{1}{2} ∫_{0}^{t_{f}}dt ∫_{Ω}dx (rn_{1}^{2}(x, t)+s (u−u_{0})^{2} (x, t))
+ ∫_{Ω}dx (q_{1}n_{1}^{2}(x, t_{f}) + q_{3}c^{2}(x, t_{f})) .  (7)

The goal is to minimize this functional with respect to the drug input rate relative to some threshold rate u_{0} and the terminal costs at t_{f}, i.e., min_{u} [J(u)]. Note that here r > 0 is the tumor burden cost coefficient and s > 0 is the drug delivery cost coefficient, while q_{1} > 0 and q_{3} > 0 are the corresponding final costs. We are trying to minimize the density of tumor cells and the drug delivery quadratic control term (u(x, t) − u_{0}(x, t))^{2}. We could have chosen a linear control which would have been less realistic, and also would give rise to problems like singular control. In addition no assumption is made about the control constraints, even though there might be physical restriction on the amount of drugs that can be administered.

IV. Vector Form

For the sake of brevity we put the mathematical model in vector form with vectors in boldface.

A. Governing equations

The vector state satisfies the PDE:

\frac{∂Y}{∂t} = D∇_{x}^{2}[Y] + Af(Y) + B(Y, t)Y + U,  (8)

where

A = [\begin{matrix}a_{1} & 0 & 0 \\ 0 & a_{2} & 0 \\ 0 & 0 & a_{3} \end{matrix}] = [a_{i}δ_{i,j}]_{3×1}, f(Y) = [\begin{matrix}f_{1}(n_{1}) \\ f_{2}(n_{2}) \\ f_{3}(c) \end{matrix}] ,

B(Y, t) = − [\begin{matrix}α_{1,2}n_{2} + κ_{1,3}c & 0 & 0 \\ 0 & α_{2,1}n_{1} + κ_{2,3}c & 0 \\ 0 & 0 & 0 \end{matrix}] ,

U(x, t) = [\begin{matrix}0 \\ 0 \\ u(x, t) \end{matrix}] = U_{3}(x, t)e_{3},  (9)

where e_{i} is the ith unit vector.

B. Objective Functional

The objective in vector form is

J[Y, U] = \frac{1}{2}∫_{0}^{t_{f}}dt ∫_{Ω}dx (Y^{⊤}RY+(U−U_{0})^{⊤}S (U−U_{0}))
+\frac{1}{2} ∫_{Ω}dx (Y^{⊤}QY)|_{t=t_{f}} ,  (10)

where R = re_{1}e_{1}^{⊤}, S = se_{3}e_{3}^{⊤}, Q = q_{1}e_{1}e_{1}^{⊤} + q_{3}e_{3}e_{3}^{⊤} and U_{0} = u_{0}(x, t)e_{3}.

V. Defining the Pseudo-Hamiltonian

We have three vectors for the Lagrange multipliers, two of which are functions of space and time and one is independent of time, needed to include the optimization constraints in the extended objective for the state PDE (8), the boundary condition (6) and the initial condition (5),

ξ(x, t) = [\begin{matrix}ξ_{1} \\ ξ_{2} \\ ξ_{3} \end{matrix}] , η(x, t) = [\begin{matrix}η_{1} \\ η_{2} \\ η_{3} \end{matrix}] , χ(x) = [\begin{matrix}χ_{1} \\ χ_{2} \\ χ_{3} \end{matrix}] , (11)

i.e., ξ_{i} = ξ_{i}(x, t), η_{i} = η_{i}(x, t) and χ_{i} = χ_{i}(x), for i = 1 : 3. Letting Z = (Y, U, ξ, η, χ), define the pseudo-Hamiltonian as,

H(Z) = \frac{1}{2}∫_{0}^{t_{f}}dt ∫_{Ω}dx (Y^{⊤}RY+(U−U_{0})^{⊤}S (U−U_{0}))

+\frac{1}{2} ∫_{Ω}dx (Y^{⊤}QY)|_{t=t_{f}}

+∫_{0}^{t_{f}}dt ∫_{Ω} dx ξ^{⊤} (\frac{∂Y}{∂t} −D∇_{x}^{2}[Y]−Af(Y)

−B(Y, t)Y − U)

+∫_{0}^{t_{f}}dt ∫_{∂Ω} dΓ η^{⊤}(−D (N̂·∇_{x}) [Y])

+ ∫_{Ω} dx χ^{⊤}(Y|_{t=0}− Y_{0}) .  (12)

VI. Optimal Control Variational Formulation

The calculus of variations is used to find differential equation of optimal control for the control, state and the co-state (adjoint or Lagrange multiplier) by seeking the functional critical point necessary conditions for the first variation [3], [8] of the pseudo-Hamiltonian H(Z).

A. Pseudo-Hamiltonian First Variation

Let the extended state vector be perturbed about the optimal trajectory Z^{∗}, so that Z = Z^{∗} + δZ, where δZ is the perturbation. Next expand the pseudo-Hamiltonian

H(Z^{∗} + δZ) = H(Z^{∗}) + δH(Z^{∗}, δZ) + O((δZ)^{2}).

Neglecting the quadratic order terms, including the 2nd variation of H, the first variation is given by functional terms linear in δZ using (12),

δH(Z^{∗}, δZ) = ∫_{0}^{t_{f}}dt ∫_{Ω}dx ((Y^{∗})^{⊤}RδY+(U^{∗}−U_{0})^{⊤}SδU)
+ ∫_{Ω}dx ((Y^{∗})^{⊤}QδY)|_{t=t_{f}}

+∫_{0}^{t_{f}}dt ∫_{Ω}dx ((ξ^{∗})^{⊤}(δY_{t}−D∇_{x}^{2}[δY]

−A(δY·∇_{Y})[f](Y^{∗}) − B(Y^{∗}, t)δY

−(δY·∇_{Y})[B](Y^{∗}, t)Y^{∗} − δU)
+δξ^{⊤}(Y_{t}^{∗}−D∇_{x}^{2}[Y^{∗}]−Af(Y^{∗})
−B(Y^{∗}, t)Y^{∗} − U^{∗}))  (13)
−∫_{0}^{t_{f}} dt∫_{∂Ω} dΓ((η^{∗})^{⊤}D (N̂·∇_{x}) [δY]
+δη^{⊤}D (N̂·∇_{x}) [Y^{∗}])

+ ∫_{Ω}dx ((χ^{∗})^{⊤}δY+δχ^{⊤}(Y|_{t=0}− Y_{0})) .

Before the critical conditions for first variation in (13) can by used to obtain the extended state equations, the functional dependence of the higher derivatives in time and state of the extended state perturbations must be eliminated on lower order terms by one or two integrations by parts, i.e., by one,

∫_{0}^{t_{f}}dt(ξ^{∗})^{⊤}δY_{t} = (ξ^{∗})^{⊤}δY|_{0}^{t_{f}} −∫_{0}^{t_{f}}dtδY^{⊤}ξ_{t}^{∗}

and by two using the Green’s formula [4],

∫_{Ω}dx(ξ^{∗})^{⊤}D∇_{x}^{2}[δY] = ∫_{Ω}dxδY^{⊤}∇_{x}^{2}[Dξ^{∗}]

+ ∫_{∂Ω}dΓ((N̂·∇_{x})[δY^{⊤}]Dξ^{∗}

−δY^{⊤}(N̂·∇_{x})[Dξ^{∗}]).

Merging these identities with (13), rearranging inner products and collecting terms the extended state equations yields the following intermediate form:

δH(Z^{∗}, δZ) = ∫_{0}^{t_{f}}dt ∫_{Ω}dx δY^{⊤}(RY^{∗}−ξ_{t}^{∗}−∇_{x}^{2}[Dξ^{∗}]

−∇_{Y} [f](Y^{∗})Aξ^{∗}−B(Y^{∗}, t)ξ^{∗} −∇_{Y}[B](Y^{∗}, t)):(ξ^{∗}(Y^{∗})^{⊤})

+∫_{0}^{t_{f}}dt ∫_{Ω}dx δU^{⊤}(S (U^{∗}−U_{0})−ξ^{∗})

+∫_{0}^{t_{f}}dt ∫_{Ω}dx δξ^{⊤}(Y_{t}^{∗}−D∇_{x}^{2}[Y^{∗}]

−Af(Y^{∗})−B(Y^{∗}, t)Y^{∗} − U^{∗})

−∫_{0}^{t_{f}}dt ∫_{∂Ω}dΓ δη^{⊤}D (N̂·∇_{x}) [Y^{∗}]

+∫_{0}^{t_{f}}dt ∫_{∂Ω}dΓ δY^{⊤}(N̂·∇_{x}) [Dξ^{∗}]

−∫_{0}^{t_{f}}dt ∫_{∂Ω}dΓ(N̂·∇_{x})[δY^{⊤}]D(η^{∗}+ξ^{∗})

+ ∫_{Ω}dx (δχ^{⊤} (Y^{∗}−Y_{0}(x)))|_{t=0}

+ ∫_{Ω}dx (δY^{⊤}(χ^{∗}−ξ^{∗}))|_{t=0}

+ ∫_{Ω}dx (δY^{⊤}(ξ^{∗}+QY))|_{t=t_{f}} ,

where A: B denotes the trace of the matrix AB or the double-dot product.

B. State Equations

The optimal state equation is recovered by setting the coefficient of (δξ^{∗})^{⊤} to zero:

\frac{∂Y^{∗}}{∂t} = D∇_{x}^{2}[Y^{∗}] + Af(Y^{∗}) + B(Y^{∗}, t)Y^{∗} + U^{∗} (14)

on Ω × (0, t_{f}], with boundary conditions on ∂Ω × [0, t_{f}] from the coefficient of (δη^{∗})^{⊤}, i.e.,

−D(N̂·∇_{x})[Y^{∗}](x, t) = 0,  (15)

for (x, t) ∈ ∂Ω × [0, t_{f}] and with initial conditions on the interior Ω from the coefficient of (δχ^{∗})^{⊤}, i.e.,

Y^{∗}(x, 0) = Y_{0}(x)  (16)

for x ∈ Ω. Due to the presence of the functions f(Y) and B(Y, t)Y the forward PDE (14) will be nonlinear.

C. Degenerate Regular Optimal Control

Since the control has been defined in (9) has been defined as only having on component, only the coefficient of δU_{3} is set to zero giving the regular, but degenerate, control as

U_{3}^{∗}(x, t) = U_{0,3}(x, t) + \frac{1}{s}ξ_{3}^{∗}(x, t),  (17)

on Ω × [0, t_{f}], provided s ≠ 0. Note that the control law with δU_{1} = 0 = δU_{2} only requires solving for the 3rd component of the first co-state vector ξ^{∗}(x, t).

D. Co-State Equations

Upon setting the functional coefficient of (δY^{∗})^{⊤} to zero yields the primary co-state backward PDE:

0 = \frac{∂ξ^{∗}}{∂t} + ∇_{x}^{2}[Dξ^{∗}] + ∇_{Y} [f](Y^{∗})Aξ^{∗}
+B(Y^{∗}, t)ξ^{∗} + ∇_{Y} [B](Y^{∗}, t)):(ξ^{∗}(Y^{∗})^{⊤})
−RY^{∗},  (18)

for (x, t) ∈ Ω × [0, t_{f}). This PDE (18) is unidirectionally coupled to the state PDE (14), except that only the 3rd component ξ_{3}^{∗}(x, t) is needed for the regular optimal control input U_{3}^{∗}(x, t) from (17). The boundary condition follows from setting the functional coefficient of δY(x, t) for x on Γ=∂Ω to zero, so

(N̂·∇_{x})[Dξ^{∗}](x, t) = 0, (x, t) ∈ ∂Ω × [0, t_{f}) (19)

and the final condition for this backward PDE follows from forcing the coefficient of δY(x, t_{f}) to be zero on Ω,

ξ^{∗}(x, t_{f}) = −QY(x, t_{f}).  (20)

The two other co-state vectors should not be needed, but satisfy rather simple equations. The 2nd co-state vector equation follows as the zero coefficient of (N̂·∇_{x})[δY^{⊤}] on the state boundary Γ=∂Ω,

η^{∗}(x, t) = −ξ^{∗}(x, t), (x, t) ∈ ∂Ω × [0, t_{f}].

The 3rd co-state vector equation follows as the zero coefficient of state initial condition δY(x, 0),

χ^{∗}(x) = ξ^{∗}(x, 0), x ∈ Ω.

VII. Iterative Approximations for Coupled System

We need to solve the system of equations developed in the previous section, namely the state equations (14) using the regular optimal control (17) and co-state equations (18), with the understanding that the state equations are forward equations while the co-state equations are backward equations. The method a double shot method, since model has two vector-valued PDEs, so the method consists of one forward shot with (14) and and one backward shot with (18).

The following steps can be used to get approximate solutions:

1) An initial guess for the first forward-backward shot iteration is made for the control U_{3}(x, t) = U_{3}^{(1)}(x, t) in (17). Substituting it into the state forward PDE (14) solving for Y^{(1)}(x, t), using the boundary condition (15) and initial condition (16).

2) Next the approximate final condition ξ^{(1)}(x, t_{f}) = −QY^{(1)}(x, t_{f}) (20) is used to start the backward co-state PDE (18) approximation using boundary conditions (19) .

3) Once ξ^{(1)}(x, t) is determined the regular optimal control equation is used to determine the updated value of the control,

U_{3}^{(2)}(x, t) = U_{0,3}(x, t) + ξ_{3}^{(1)}(x, t)/s.

4) This process is repeated for ℓ = 2 : L double shot iterations until a convergence criterion for sufficiently large L is reached, e.g., the relative criterion for the control,

||U_{3}^{(ℓ)}(x, t)−U_{3}^{(ℓ−1)}(x, t)||<tol_{u}||U_{3}^{(ℓ−1)}(x, t)|| ,

and, say,

||Y_{1}^{(ℓ)}(x, t)−Y_{1}^{(ℓ−1)}(x, t)||<tol_{y}||Y_{1}^{(ℓ−1)}(x, t)|| ,

for ℓ = 2:L until satisfied, provided ||U_{3}^{(ℓ−1)}(x, t)|| ≠ 0 and ||Y_{1}^{(ℓ−1)}(x, t)|| ≠ 0, where tol_{u} > 0 and tol_{y} > 0 are some prescribed tolerances.

VIII. Double Shot Forward-Backward Computational Method

In reality the problem is highly non-linear as are many problems in biology and we need numerical approximations of the solution. The main problem here is the fact that we have a forward state equation and a backward co-state equation. What we really do here is modify shooting methods [7] for initial-final-boundary value problems, where the starting aim is replaced by an estimate of the full control law U_{3}^{(ℓ)}(x, t) for the forward integration of the state PDE (14) whose final approximation Y^{(ℓ)}(x, t_{f}) serves as the backward aim (20) for the backward integration of the costate PDE (18) producing an approximation ξ^{(ℓ)}(x, t) whose third component is used to update (17) the control law U_{3}^{(ℓ+1)}(x, t).

Before we actually go into the details of the shooting method, we will briefly discuss the Crank-Nicolson, predictor-corrector central finite difference method used to solve the nonlinear PDEs. The space is discretized as follows,

x → x_{j} = [x_{j_{i},1} + (j_{i} − 1) · Δx_{i}]_{3×1}.

Here Δx_{i} is the mesh size for state i and j = [j_{i}]_{3×1} where, j_{i} = 1:M_{i} nodes per node for states i = 1:3. For the forward state equation we have the forward time discretization, t → t_{k} = kΔt, for k = 0:K time steps where Δt is the forward time step size, t_{0} = 0 and t_{K} = t_{f}. Now we consider the vector state/co-state PDE system in general form,

Y_{t}^{∗} = F(x, t, Y^{∗}(x, t), U^{∗}(x, t)),
0 = ξ_{t}^{∗} + G(x, t, ξ^{∗}(x, t), Y^{∗}(x, t)).

In this equation, using the ordinate discretization Y(x_{j}, t_{k}) ≃ Y_{j,k} and ξ(x_{j}, t_{k}) ≃ ξ_{j,k}, the relevant derivatives discretization are

(Y_{t}^{∗})_{j,k+1/2}≃\frac{Y_{j,k+1/2+1/2} − Y_{j,k+1/2−1/2}}{2(Δt/2)} ,
(ξ_{t}^{∗})_{j,k−1/2}≃\frac{ξ_{j,k−1/2+1/2} − ξ_{j,k−1/2−1/2}}{2(Δt/2)}

and

(∇_{x}^{2}[Y^{∗}])_{j,k+1/2}≃\frac{Y_{j+e_{i},k+1/2}−2Y_{j,k+1/2}+Y_{j−e_{i},k+1/2}}{(Δx_{i})^{2}} ,
(∇_{x}^{2}[ξ^{∗}])_{j,k−1/2}≃\frac{ξ_{j+e_{i},k−1/2}−2ξ_{j,k−1/2}+ξ_{j−e_{i},k−1/2}}{(Δx_{i})^{2}} ,

Consequently the forward and backward numerical schemes are given by

Y_{j,k+1}^{(ℓ)} = Y_{j,k}^{(ℓ)} + ΔtF_{j,k+ \frac{1}{2}}^{(ℓ)} , for k = 0:K−1,
ξ_{j,k−1}^{(ℓ)} = ξ_{j,k}^{(ℓ)} + ΔtG_{j,k− \frac{1}{2}}^{(ℓ)} , for k =:−1:1,

respectively, for ℓ = 1:L. For each double shot ℓ the state starts from Y_{j,0}^{(ℓ)} = Y_{0,j,0} = Y_{0}(x_{j}, 0) using

U_{3,j,k}^{(ℓ)} = U_{0,3,j,k} + ξ_{3,j,k}^{(ℓ−1)}/s, for k = 0:K−1,

except when ℓ = 1 and the initial guess U_{3,j,k}^{(1)}, e.g., U_{0,3,j,k}, is used. For each updated forward state shot is completed, then the backward co-state shot starts from ξ_{j,K}^{(ℓ)} = QY_{j,K}^{(ℓ)} using the whole state set Y_{j,k}^{(ℓ)} for k = 0 : K. The Crank-Nicolson temporal mid-point F_{j,k+ \frac{1}{2}}^{(ℓ)} state function is approximated by average with

Y_{j,k+ \frac{1}{2}}^{(ℓ)} = \frac{1}{2} (Y_{j,k+1}^{(ℓ)} + Y_{j,k}^{(ℓ)}) ,

which can be used to construct finite differences for the derivatives, with a similar form for U_{j,k+ \frac{1}{2}}^{(ℓ)} , for ℓ = 0:L.

Similarly, the backward temporal mid-point G_{j,k− \frac{1}{2}}^{(ℓ)} co-state function is approximated with the average,

ξ_{j,k− \frac{1}{2}}^{(ℓ)} = \frac{1}{2} (ξ_{j,k}^{(ℓ)} + ξ_{j,k−1}^{(ℓ)}) ,

and is used for the derivatives as well.

The no flux boundary conditions for both the state and co-state present some extra complexity, since the central differences of Crank-Nicolson are not suitable at the boundary if it is necessary to avoid using artificial external points. External points can be avoided by judicious use of forward and backward differences of second order, matching the accuracy of the Crank-Nicolson central differences. In the simplest case of rectangular grids, the discretized no flux boundary conditions (15,19) with second order accuracy are

0 = − ((N̂·∇_{x})[Y^{∗}])_{j,k}^{(ℓ)}

≃ −\frac{3Y_{j,k}^{(ℓ)} − 4Y_{j−N,k}^{(ℓ)} + Y_{j−2N,k}^{(ℓ)}}{2|N·Δx|} ,
0 = ((N̂·∇_{x})[(ξ)^{∗}])_{j,k}^{(ℓ)}

≃ \frac{3(ξ)_{j,k}^{(ℓ)} − 4(ξ)_{j−N,k}^{(ℓ)} + (ξ)_{j−2N,k}^{(ℓ)}}{2|N·Δx|} ,

respectively, where N = N̂_{j,k}, Δx = [Δx_{i}]_{3×1} > 0, D is not needed, and, e.g., Y_{j−N,k}^{(ℓ)}=Y^{(ℓ)}(x_{j} − |N·Δx|N, t_{k}). For non-rectangular domains, interpolation would be needed to convert evaluations to defined spatial nodes or else domain compatible grids should be used, e.g., for circular or spherical grid boundaries, N = e_{r}, where here r is the radius and N · Δx = Δr.

During each ℓth double shot, a prediction and corrections of the state and co-state are used to account for the usual nonlinearities in the biological models [6], stopping when the changes are sufficiently small.

The above mentioned method is a sequential double shot method since one shot is used to get Y_{j,k}^{(ℓ)} and a subsequent shot it used to get ξ_{j,k}^{(ℓ)}. Alternately, a parallel two shot method could be used to get an approximate solution by integrating for both Y_{j,k}^{(ℓ)} and ξ_{j,k}^{(ℓ)} in the forward direction using a guess initial condition for ξ_{j,0}^{(ℓ)} at t_{0} = 0, with several genuine shooting method shots until some ℓ^{∗} shot where ||ξ_{j,K}^{(ℓ^{∗})} +QY_{j,K}^{(ℓ^{∗})}|| < tol_{ξ}, i.e., the final co-state value is small enough using some sufficiently small tolerance tol_{ξ} to approximate the final condition (20). Interpolation should be used to accelerate the convergence of ξ_{j,0}^{(ℓ)} to some sufficiently small value.

IX. Test Results

The double shot forward and backward algorithm outlined in the previous two sections has been tested on one space dimension, x, example with three state dimensions {Y_{1}^{∗}=N_{1}^{∗}, Y_{2}^{∗}=N_{2}^{∗}, Y_{3}^{∗}=C^{∗}}, plus the drug input control U_{3}^{∗}. The numerical parameter data come from the BCNU drug simulations for the brain of Wang et al. [14] and the brain tumor modeling of Swanson [12] and Murray [11], with some difficult to find parameters from Mansuri [9] or from reasonable estimates from other areas. For example, diffusion diagonal vector is D = [4.2e-3, 1.e-15, 2.2e-1] cm^{2} per day (normal tissue diffusion is assumed to be insignificant), the quadratic cost coefficients are r = 0.1 = s = q_{1} = q_{3}, the net growth coefficients are a = [1.2e-2, 8.6e-7, 1.1e+1] per day, the carrying capacities for tumor k_{1} and normal k_{2} tissues are scaled to one for the normal value and the interaction coefficients {α_{1,2}, α_{2,1}, κ_{1,3}, κ_{2,3}} are all given the arbitrary value 1.0e-4. The initial states are given to be uniformly one for the normal tissue, while the tumor density was assumed to be a spatial Gaussian with spread 2.0e-2 about a mean of 0.0 with a weight of 1.e-3. The initial drug concentration has a Gaussian spread of 1.0e-2 about a mean 0.0 with weight 1.e-4, while the threshold drug control u_{0}(x, t) is similarly distributed, but with weight 1.e-5.

The results are given only for the tumor density N_{1}^{∗} = N_{1}^{∗}(x, t) in Figure 1 on the symmetric space interval x ∈ [−1.0, +1.0] in centimeters over a t_{f} = 5 day treatment. For this simple one space dimension test example, we see that the optimal distribution of the tumor using an optimal distribution of the drug delivery results in the 16.2% reduction of the tumor density over this simulated five day drug treatment trial.

Fig. 1. Tumor density N_{1}^{∗}(x, t) versus the one-dimensional spatial coordinate x with time t at the rounded parameter values {0, t_{q1} = 0.25t_{f} , t_{mid} = 0.5t_{f} , t_{q3} = 0.75t_{f} , t_{f} }, where t_{f} = 5 days. The targeted tumor density rapidly decays in this simulated 5 day trial.

X. Conclusion

The main interest of this paper was to provide the necessary foundation to study the mechanism of drug delivery to the brain. We have set up a fairly realistic distributed parameters model which takes into account the spatial dependence of the state variables. The main focus of the paper was to develop an algorithm to determine the optimal drug delivery to brain tumors using an optimal distribution of the drug about the original tumor site. This paper leaves room for a lot of new directions of work.

A. Future Directions One such direction would be running a simulation to

implement the algorithm, using supercomputing tools. Of course this would require more realistic medical data. Another important aspect that can be examined is the effect of stochasticity, most notably the Gaussian and Poisson type of noise. The physical basis for such stochasticity would be the phenomenon of metastasis, which gives rise to additional tumor sites and also the side effects produced by the drug. We hope to examine all these aspects of the problem in future works.

XI. Acknowledgments

The authors gratefully acknowledge, Dr. Herbert H. Engelhard of Neurological Surgery at the University of Illinois Medical Center for suggesting the problem and Professor William D. O’Neill of Bioengineering at the University of Illinois at Chicago for putting us in contact.

References

[1] H.H. Engelhard, “Brain Tumors and the Blood-Brain Barrier,” Neuro-Oncology, The Essentials, pp. 49-53, Thieme Medical Publishers, Inc., New York 2000.

[2] R.A. Gatenby and E.T. Gawlinski, “A Reaction-Diffusion Model of Cancer Invasion,” Cancer Research, vol. 56, pp. 5745-5753, December 15, 1996.

[3] M.D. Gunzburger, Perspectives in Flow Control and Optimization, SIAM Series in Design and Control, Philadelphia, 2003.

[4] R. Haberman, Elementary Applied Partial Differential Equations with Fourier Series and Boundary Value Problems, Prentice-Hall, Englewood Cliffs, NJ, 1983.

[5] F.B.Hanson and K.Naimipour, “Convergence of Numerical Method for Multistate Stochastic Dynamic Progamming,” Proceedings of International Federation of Automatic Control 12th World Congress, vol. 9, pp.501-504, 1993.

[6] F.B.Hanson, “Techniques in Computational Stochastic Dynamic Programming,” Control and Dynamic Systems: Advances in Theory and Applications, vol. 76, pp.103-162, Academic Press, New York, April 1996.

[7] H. B. Keller, Numerical Methods for Two-Point Boundary Value Problems, Gin Blaisdell, Waltham, MA, 1968.

[8] D. E. Kirk, Optimal Control Theory: An Introduction, Prentice-Hall, Englewood Cliffs, NJ, 1970.

[9] J. Mansuri, The Modelling of Tumor Growth Using Reaction-Diffusion Equations, MSc Thesis, Oxford University, 2002.

[10] J.D.Murray, Mathematical Biology I: An Introduction, Springer, New York, 2002.

[11] J.D.Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, Springer, New York, 2003.

[12] K.R. Swanson, Mathematical Modeling of the Growth and Control of Tumors, PhD Thesis, University of Washington, 1999.

[13] C-H. Wang and J. Li, “Three dimensional simulation of IgG delivery to tumors,” Chemical Engineering Science, vol. 53, pp. 3579-3600, 1998.

[14] C-H. Wang, J. Li, C.S. Teo and T. Lee, “The delivery of BCNU to brain tumors,” Journal of Controlled Release, vol. 61, pp. 21-41, 1999.

[15] J. J. Westman, B. R. Fabijonas, D. L. Kern, and F. B. Hanson, “Compartmental Model for Cancer Evolution: Chemotherapy and Drug Resistance,” Mathematical BioSciences, pp. 1-20, May 2002, under revision.

[16] D.E.Woodward, J.Cook, P.Tracqui, G.C.Cruywagen, J. D Murray, E.C.Alvord, Jr., “A mathematical model of glioma growth: the effect of extent of surgical resection,” Cell Proliferation, vol. 29, pp. 269-288, 1996.