Information Technology Journal

Year: 2010 | Volume: 9 | Issue: 8 | Page No.: 1647-1652
DOI: 10.3923/itj.2010.1647.1652

A Computer Algorithm for Optimizing to Extract Effective Diffusion Coefficients of Drug Delivery from Cylinders

Youyun Li, Ziquan Xiang, Xuyu Xiang and Song Wang

Abstract: In this study, we propose a computer optimization algorithm for estimating effective diffusion coefficients of drug delivery from a cylinder to an external cylindrical finite volume. We first write the diffusion equation in the polar-coordinate form and then a finite difference scheme for the diffusion equation is developed for solving the equation. The diffusion coefficient extraction is formulated as a least squares problem. To solve the lest-squares problem defining the unknown diffusivity, a computer algorithm of Gauss-Newton type is proposed. Numerical results are presented to validate the numerical methods proposed.

Keywords: Gauss-Newton scheme, Effective diffusion coefficient and algorithm for optimization

INTRODUCTION

Diffusion and convection-diffusion processes appear in many areas such as geo-physics, engineering, biomedical science (Mwellott et al., 2001; Fu et al., 1976; Grassi and Grassi, 2005; Siepmann et al., 1998; Crawford et al., 2002; Hicks et al., 2003). In many cases, diffusion coefficients, or diffusivity, are unknown and need to be identified using experimentally or exploratory observed data. While a diffusion process can be governed by a function of space, time and concentration of substance, in practice, we normally find a constant diffusion coefficient to approximate the process. The design of controlled drug delivery devices has attracted much of attention for that which the effective diffusivity of a device is critical to its functionality and performance (Lou et al., 2005; Price Jr. et al., 1997; Asaoka and Hirano, 2003; Hukka, 1999; Kohne et al., 2002). Although the diffusivity of a drug delivery system is determined mainly by the porosity and some other properties of the materials (Lou et al., 2005; Price Jr. et al., 1997; Asaoka and Hirano, 2003), when these properties are known, how to extract the effective diffusivity of the system becomes a major concern. There are various existing techniques for the identification of effective diffusivity. These techniques are based on either empirical or semi-empirical models from drug delivery mechanism or on analytical solutions of the diffusion equation in 2D or in the special cases (Price Jr. et al., 1997; Asaoka and Hirano, 2003; Kohne et al., 2002). In practical application, devices are always three-dimensional. It is difficult to extract the diffusion parameters if only depending on the empirical or semi-empirical models. However, for the analytical solution, the cases only are limited in 2D or special devices (Wang and Lou, 2007; Lou et al., 2004). Therefore, in order to better analyze 3D cases, it is necessary to establish new numerical methods to extract the diffusion coefficients from the diffusion and convection-diffusion processes. In this study, we shall propose a numerical method based on a finite difference scheme to estimate the diffusion coefficient from the 3D drug delivery system. Finally, numerical results are shown to illustrate the convergence and usefulness of the optimization method.

THE PROBLEM

Consider a device of cylinder with radius r₁ and height h₁ loaded an amount M⁰ of drug. This device is placed in a cylindrical container of radius r₂ and height h₂ filled with liquid. The configuration is shown in Fig. 1. We let the regions of the device and the container be denoted respectively by Ω₁ and Ω. The diffusion process of this problem with a constant coefficient is governed by the following equation in Cartesian coordinates:

where, D is a constant and C(x, y, z, t) is the unknown concentration.

For the initial condition H(x, y, z), we assume that at t = 0, the concentration is uniform in the device and zero in liquid, i.e.:

Therefore, the process to determine the diffusion coefficient D is equivalent to solving the following optimization problem.

Problem 1: Find D to satisfy:

where, M₁⁰, M₂⁰, ..., M_e⁰ are given experimental data and M₁, M₂, ..., M_e are computed using the following expression:

with C(x, y, z, t) being the solution to Eq. 1:

To compute C, we shall compute Problem 1 in the polar coordinates. Therefore, Problem 1 can be transformed into the following Problem 2.

Problem 2: Find D to satisfy:

where, M₁, M₂, ..., M_e are computed as the following equations:

and C(r, θ, z, t) is governed by the following polar coordinate diffusion equation:

with the initial condition:

DISCRETIZATION

To determine the diffusion coefficients D, it is necessary to solve the partial differential Eq. 7. Here, we propose a finite difference scheme for the discretization of Eq. 7.

Let the intervals (0, 2π), (0, r₂) and (0, h) be divided uniformly into P, Q and R sub-intervals, respectively with the respective interval lengths Δr, θ and Δz, where, P, Q and R are given positive integers. This defines a mesh for the container region with mesh nodes (θ_i, r_j, z_k) for i = 1 ..., P, j = 1, ..., Q and k = 1, ..., R. For a given time step length Δt and a given positive integer T, we let l = 1, 2, ..., T for t_l = (l-1)Δt. Using this partition in space and time, we define the following approximations for the derivatives appearing in Eq. 7:

Using Eq. 9-13 and 7 can be discretized into the following difference equation system:

for all admissible C^l_i,j,k, where, D₁, D₂, ..., D₇ are defined by:

This defines a linear system for the unknowns C^l_{i, j, k} for all admissible i, j, k, l.

We comments though a uniform partition is used in the above discussion for brevity of notation, it is obvious that the discretization scheme is also true for non-uniform partitions.

NUMERICAL METHODS FOR SOLVING THE LEAST-SQUARES PROBLEM

In this section, we will present some algorithms for Problem 2. Let:

where, M = (M₁ (D), M₂ (D), ..., M_e (D))^T and M* = (M⁰₁, M⁰₂, ..., M⁰_e)^T.

We now consider the numerical solution of Problem 2. Starting from an initial guess D⁰, Problem 2 can be solved iteratively. At each step an increment δDⁱ is calculated such that:

is minimized with respect to δDⁱ, Dⁱ and δDⁱ are the ith approximation and ith increment of D, respectively. The iterative procedure continues until the relative error:

is smaller than a given small positive constant.

The gauss-newton method: To calculate the increment δDⁱ at each step, in this study, based on the idea given by Lee et al. (1999), one Gauss-Newton method is established.

Taylor's formula for vector valued functions gives:

where,

and G denotes the second derivative vector of Mⁱ evaluated at Dⁱ+ρδDⁱ with 0≤ρ≤, in this study, we set ρ = 1. Omitting the second order terms in Eq. 17, we have:

When δDⁱ is small, E(Dⁱ+δDⁱ) can be approximated by:

This is a quadratic form in δDⁱ and the minimum point δDⁱ* of this quadratic function satisfies:

which leads to:

The solution to Eq. 20 defines the ith search direction called the Gauss-Newton direction. Solving Eq. 20 gives:

The new approximation to the diffusion coefficients D is then defined as:

Evaluation of partial derivatives in J_i: In order to obtain the ith search direction δDⁱ*, from Eq. 21, it is necessary to compute all the partial derivatives in J_i for Dⁱ. In what follow, we shall present an algorithm to compute J_i.

From Eq. 6, we get the following equation:

In order to obtain the derivative using Eq. 23, we firstly compute the derivatives For any computed diffusion coefficient Dⁱ, we let:

where, δ is a small constant. Based on the diffusion Eq. 1 and the ith step diffusion coefficients Dⁱ and DDⁱ, we can compute the concentration C_i of Dⁱ and the concentration C_ii of DDⁱ. Using the following difference formula:

the derivative can be approximated. Therefore, from Eq. 23, the derivative can be obtained as follows:

We comment that the above partial derivatives are computed by the finite difference method shown as above. If the discrete value of the partial derivatives are obtained, the discrete value of the partial derivatives can be approximated by the Eq. 23.

THE LEAST-SQUARES COMPUTER ALGORITHM

The following algorithm is based on the numerical methods presented in the previous sections.

Let M_i, i = 1, 2, 3, .., e be a set of the mass points by the numerical equation and M_i* be a set of the experimentally measured mass at t_i for each i = 1, 2, 3, .., e. The diffusion coefficient D can be determined by the following least squares (LSQ) algorithm.

ALGORITHM LSQ

NUMERICAL EXAMPLES

To verify the usefulness of the numerical methods, numerical experiments were performed. In the numerical experiments, non-uniform partitions are used.

Only testing the validity of the computer algorithm, for the brief, we only introduce the mathematical and computer model, not considering the material. The test problem done by the school of mathematics and statistics in UWA in 2007 is a cylindrical device placed in a cylindrical container shown in Fig. 1 with their sizes given in Table 1. Table 2 lists The experimental release data (M_t*/M_∞) at the different time points in the lab. In order to reduce the computing time, we first solve this optimization problem using the initial starting point D₀ = 1x10^-1 using the mesh parameters Δt = 100s, Δr = 0.002, Δr = 0.002 cm, Δz = 0.04 cm and Δθ = 0.8 cm. Table 3 gives the computed data (M_t/M_∞) from the first 6 iterations and the last 4 iterations by the computer algorithm. We then reduce the mesh sizes to Δt = 50s, Δr = 0.001 cm, Δz = 0.02 cm and Δθ = 0.8 cm and use the results from the 6th iteration, D₀ = 2.25276x10^-6 and the initial guess to continue our computation for 4 more iterations.

Table 4 shows iterated step distances and the errors for the drug releasing delivery system by the optimal algorithm. In the numerical experiments, we set the error in each step as:

Table 4 shows the step distances of the computer searching algorithms. Table 5 gives the total error between the lab data and the computed data by the computer algorithm. Figure 2 plots the computed diffusion coefficients. Figure 3 plots the release profile. From Table 4, 5, Fig. 2 and 3, the optimal algorithm based on FD is convergence. From the last lines in Table 5, the least-squares error is small. Figure 4 plots the fitted curves which also indicated that the computed results (M_t/M_∞) are very close to the experiment data (M_t*/M_∞). Therefore, from the analysis of the computed data, the optimization algorithm based on FD is valid.

CONCLUSION

In this study, we developed some mathematical optimal numerical methods based on the finite difference method for estimating effective diffusiveness of a drug from a delivery device of tube geometry in three dimensions to an external finite volume. Numerical experiments were performed and the numerical results show the usefulness of the methods developed.

ACKNOWLEDGMENTS

This study is supported by a project supported by Scientific Research Foundation for Returned Scholars, Ministry of Education of China and by the National Natural Science Foundation of China (50778026), an Australian Research Council Discovery Grant (DP0557148).

REFERENCES