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
How to cite this article
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.
Consider a device of cylinder with radius r1 and height h1 loaded an amount M0 of drug. This device is placed in a cylindrical container of radius r2 and height h2 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:
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, M10, M20, ..., Me0 are given experimental data and M1, M2, ..., Me 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, M1, M2, ..., Me are computed as the following equations:
and C(r, θ, z, t) is governed by the following polar coordinate diffusion equation:
with the initial condition:
Let the intervals (0, 2π), (0, r2) 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, rj, zk) 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 tl = (l-1)Δt. Using this partition in space and time, we define the following approximations for the derivatives appearing in Eq. 7:
for all admissible Cli,j,k, where, D1, D2, ..., D7 are defined by:
This defines a linear system for the unknowns Cli, 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 = (M1 (D), M2 (D), ..., Me (D))T and M* = (M01, M02, ..., M0e)T.
We now consider the numerical solution of Problem 2. Starting from an initial guess D0, Problem 2 can be solved iteratively. At each step an increment δDi is calculated such that:
is minimized with respect to δDi, Di and δDi 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 δDi 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:
and G denotes the second derivative vector of Mi evaluated at Di+ρδDi with 0≤ρ≤, in this study, we set ρ = 1. Omitting the second order terms in Eq. 17, we have:
M = Mi+JiδDi
When δDi is small, E(Di+δDi) can be approximated by:
This is a quadratic form in δDi and the minimum point δDi* of this quadratic function satisfies:
ΔE(Di+δDi) = 0
which leads to:
The new approximation to the diffusion coefficients D is then defined as:
Evaluation of partial derivatives in Ji: In order to obtain the ith search direction δDi*, from Eq. 21, it is necessary to compute all the partial derivatives in Ji for Di. In what follow, we shall present an algorithm to compute Ji.
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 Di, we let:
DDi = Dix(1+δ)
where, δ is a small constant. Based on the diffusion Eq. 1 and the ith step diffusion coefficients Di and DDi, we can compute the concentration Ci of Di and the concentration Cii of DDi. 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 Mi, i = 1, 2, 3, .., e be a set of the mass points by the numerical equation and Mi* be a set of the experimentally measured mass at ti for each i = 1, 2, 3, .., e. The diffusion coefficient D can be determined by the following least squares (LSQ) algorithm.
|Step 1:||Choose a positive integer N, let Eopt = 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 [tj-1, tj], go to step 3, otherwise, go to step 5|
|Step 3:||Compute the derivatives of the concentration in the time interval [tj-1, tj], go to step 4|
|Step 4:||Compute the mass Mj 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)>Eopt, go to step 9, otherwise, go to step 7|
|Step 7:||Compute the derivative vector (∂Mi/∂D), formulate the derivative vector Ji and get the ith search direction δDi* 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|
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 (Mt*/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 D0 = 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 (Mt/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, D0 = 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 (Mt/M∞) are very close to the experiment data (Mt*/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 (Mt*/M∞) for the different time in the drug releasing delivery system|
|Table 3:||Computed releasing data (Mt/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 (Mt*/M∞) and the computed data (Mt/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|
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.
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).
- 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.
- Grassi, M. and G. Grassi, 2005. Mathematical modelling and controlled drug delivery: Matrix systems. Curr. Drug Deliv., 2: 97-116.
- 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.
- Lou, X., S. Vijayasekaran, R. sugiharti and T. Robertson, 2005. Morphological and topographic effect on calcification tendency of pHEMA hydrogels. Biomaterials, 26: 5808-5871.
- Asaoka, K. and S. Hirano, 2003. Diffusion coefficient of water through dental composite resin. Biomaterials, 24: 975-979.
- Hukka, A., 1999. The effective diffusion coefficient and mass transfer coefficient of nordic softwoods as calculated from direct drying experiments. Holzforschung, 53: 534-540.
- 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.
- Lou, X., S. Munro and S. Wang, 2004. Drug release characteristics of phase separation pHEMA sponge materials. Biomaterials, 25: 5071-5080.
- 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.