**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.

PDF Abstract XML References Citation

**Received:**May 20, 2010;

**Accepted:**July 12, 2010;

**Published:**August 21, 2010

####
**How to cite this article**

*Information Technology Journal, 9: 1647-1652.*

**DOI:**10.3923/itj.2010.1647.1652

**URL:**https://scialert.net/abstract/?doi=itj.2010.1647.1652

**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_{1} and height h_{1} loaded an amount M^{0} of drug. This device is placed in a cylindrical container of radius r_{2} and height h_{2} filled with liquid. The configuration is shown in Fig. 1. We let the regions of the device and the container be denoted respectively by Ω_{1} and Ω. The diffusion process of this problem with a constant coefficient is governed by the following equation in Cartesian coordinates:

(1) |

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.:

(2) |

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

**Problem 1:** Find D to satisfy:

(3) |

where, M_{1}^{0}, M_{2}^{0}, ..., M_{e}^{0} are given experimental data and M_{1}, M_{2}, ..., M_{e} are computed using the following expression:

(4) |

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:

(5) |

where, M_{1}, M_{2}, ..., M_{e} are computed as the following equations:

(6) |

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

(7) |

with the initial condition:

(8) |

**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_{2}) 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:

(9) |

(10) |

(11) |

(12) |

(13) |

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

(14) |

for all admissible C^{l}_{i,j,k}, where, D_{1}, D_{2}, ..., D_{7} are defined by:

(15) |

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:

(16) |

where, M = (M_{1} (D), M_{2} (D), ..., M_{e} (D))^{T} and M* = (M^{0}_{1}, M^{0}_{2}, ..., M^{0}_{e})^{T}.

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

E(D ^{i}+δD^{i}) |

is minimized with respect to δD^{i}, D^{i} and δD^{i} 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^{i} 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:

(17) |

where,

(18) |

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

M = M ^{i}+J_{i}δD^{i} |

When δD^{i} is small, E(D^{i}+δD^{i}) can be approximated by:

(19) |

This is a quadratic form in δD^{i} and the minimum point δD^{i}* of this quadratic function satisfies:

ΔE(D ^{i}+δD^{i}) = 0 |

which leads to:

(20) |

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

(21) |

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

(22) |

**Evaluation of partial derivatives in J _{i}:** In order to obtain the ith search direction δD

^{i}*, from Eq. 21, it is necessary to compute all the partial derivatives in J

_{i}for D

^{i}. In what follow, we shall present an algorithm to compute J

_{i}.

From Eq. 6, we get the following equation:

(23) |

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

DD ^{i} = D^{i}x(1+δ) |

where, δ is a small constant. Based on the diffusion Eq. 1 and the ith step diffusion coefficients D^{i} and DD^{i}, we can compute the concentration C_{i} of D^{i} and the concentration C_{ii} of DD^{i}. Using the following difference formula:

(24) |

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

(25) |

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**

Step 1: | Choose a positive integer N, let E_{opt} = cons tant, set i = 1 and j = 1 and give the intial value of D, then go to step 2 |

Step 2: | If j≤e, solve the diffusion equation by the formula Eq. 14 to obtain the concentration in the time interval [t_{j-1}, t_{j}], go to step 3, otherwise, go to step 5 |

Step 3: | Compute the derivatives of the concentration in the time interval [t_{j-1}, t_{j}], go to step 4 |

Step 4: | Compute the mass M_{j} using the diffusion coefficient D by the formula (6), set j = j+1, go to step 2 |

Step 5: | Based on Eq. 23, compute the derivatives of the diffusion mass for the diffusion coefficient, go to step 6 |

Step 6: | Compute the E(D) by Eq. 16, if E(D)>E_{opt}, go to step 9, otherwise, go to step 7 |

Step 7: | Compute the derivative vector (∂M_{i}/∂D), formulate the derivative vector J_{i} and get the ith search direction δD^{i}* by the formula (21), go to step 8 |

Step 8: | Set D = D+δD*^{i}, i = i +1 and j = 1, go to step 2 |

Step 9: | Output the optimal diffusion parameter D and stop |

**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_{0} = 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_{0} = 2.25276x10^{-6} and the initial guess to continue our computation for 4 more iterations.

Fig. 1: | The device of the tube container with small cylindrical tube |

Fig. 2: | Diffusion coefficient extracting from a drug in the cylinder by the optimal algorithm based on FD |

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.

Table 1: | Sizes of drug releasing delivery system |

Table 2: | Experimental releasing data (M_{t}*/M_{∞}) for the different time in the drug releasing delivery system |

Table 3: | Computed releasing data (M_{t}/M_{∞}) for the different time in the drug releasing delivery system from No.1-12 by the computer algorithm |

Table 4: | Iterated step distances and the errors for the drug releasing delivery system by the optimal algorithm |

Fig. 3: | Step distances in each iterated process by the optimal algorithm based on FD |

Fig. 4: | The comparison between the experiment data (M _{t}*/M_{∞}) and the computed data (M_{t}/M_{∞}) in iterated No.10 by the optimal algorithm based on FD |

Table 5: | Iterated step distances and the errors for the drug releasing delivery system by the optimal algorithm with smaller sizes |

**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**

- Fu, J.C., C. Hagemeir and D.L. Moyer, 1976. An unified mathematical model for diffusion from drugpolyner composite tablets. J. Biomed. Matter. Res., 10: 743-758.

PubMed - Grassi, M. and G. Grassi, 2005. Mathematical modelling and controlled drug delivery: Matrix systems. Curr. Drug Deliv., 2: 97-116.

Direct Link - Siepmann, J., A. Ainaoui, J.M. Vergnaud and R. Bodmeier, 1998. Calculaiton of the dimensions of drug polymer devices based on diffusion parameters. J. Pharm. Sci., 87: 827-832.

CrossRefDirect Link - Crawford, G.J., C.R. Hicks, X. Lou, S. Vijayasekaran, D. Tan, T.V. Chirila and I.J. Constable, 2002. The chirila keratoprostesis: Phase I human clinical trials. Ophthalmology, 109: 883-889.

PubMed - Hicks, C.R., G.J. Crawford, X. Lou, D.T. Tan and G.R. Snibson
*et al*., 2003. Corneal replacement using a synthetic hydrogel cornea, AlphaCor™: device, preliminary outcomes and complications. Eye, 17: 385-392.

CrossRefDirect Link - Lou, X., S. Vijayasekaran, R. sugiharti and T. Robertson, 2005. Morphological and topographic effect on calcification tendency of pHEMA hydrogels. Biomaterials, 26: 5808-5871.

PubMed - Price Jr., P.E., S. Wang and H.H. Romdhane, 1997. Extracting effective diffusion parameters from drying experiments. AICHe J., 43: 1925-1934.

CrossRefDirect Link - Asaoka, K. and S. Hirano, 2003. Diffusion coefficient of water through dental composite resin. Biomaterials, 24: 975-979.

CrossRef - Hukka, A., 1999. The effective diffusion coefficient and mass transfer coefficient of nordic softwoods as calculated from direct drying experiments. Holzforschung, 53: 534-540.

Direct Link - Kohne, J.M., H.H. Gerke and S. Kohne, 2002. Effective diffusion coefficieents of soil aggregates with surface skins. Soil Sci. Soc. Am. J., 66: 1430-1438.

Direct Link - Lou, X., S. Munro and S. Wang, 2004. Drug release characteristics of phase separation pHEMA sponge materials. Biomaterials, 25: 5071-5080.

CrossRef - Lee, W.R., S. Wang and K.L. Teo, 1999. An optimization approach to a finite dimensional parameter estimation problem in semiconductor device design. J. Comput. Phys., 156: 241-256.

Direct Link