Numerical methods for estimating effective diffusion coefficients of three-dimensional drug delivery systems

This paper presents a numerical technique in three dimensions for estimating effective 
diffusion coefficients of drug release devices in rotating and flow-through fluid systems. 
We first formulate the drug release problems as diffusion equation systems with unknown effective 
diffusion coefficients. We then 
develop a numerical technique for estimating the unknown coefficients based on a nonlinear least-squares method and a 
finite volume discretization scheme for the 3D diffusion equations. Numerical experiments have been performed 
using experimental data and the numerical results are presented to 
show that our methods give accurate diffusivity estimations 
for the test problems.


1.
Introduction. Many real world phenomena are governed by diffusion processes. It is very often that the diffusion coefficient of such a process is unknown and needs to be determined in order to understand the process(cf., for example, [1,22,6,8,2,21,18,20]). One example is the controlled drug delivery in which a drug is delivered to a location at a desired rate from a diffusional device. Drug release devices are usually made of bio-compatible polymeric materials. The diffusion coefficient of such a material is a measure of the porosity and other properties of the polymer matrix which determine the drug transport in the matrix. Therefore, the estimation of the diffusion coefficient is crucial for designing a controlled drug release device so that it can achieve a desired release rate. In theory, a diffusion coefficient can be a function of space, time, concentration and even temperature of the system [6]. However, in practice we usually seek for a constant approximation to the coefficient. This constant approximation is called the effective diffusion coefficient.
Given an experimentally observed time series of drug release profile of a device, the determination of the effective diffusion coefficient can be posed as a nonlinear optimization problem subject to a set of constraints involving diffusion equations. Methods for estimating effective diffusion coefficients can be found in the open literature which involve optimization techniques and solutions of diffusion equations. Analytical solutions to diffusion problems are available for some problems with simple geometries [3]. However, diffusion problems of practical significance can not usually be solved analytically. Empirical diffusion models have also been proposed to estimate effective diffusion coefficients, but they do not have solid mathematical backgrounds. Recently, we have proposed a full numerical method for estimating effective diffusion coefficients in 2-dimensions [19]. Since drug release problems in practice are always 3-dimensional, full numerical methods in 3D are desirable for solving realistic problems.
In this paper we develop a 3-dimensional full numerical technique for drug delivery problems in the rotating fluid system and the flow-through system. This technique is based on an optimization algorithm for the nonlinear least-squares problem and a discretization scheme for the diffusion equations. For discussion simplicity, we use cylindrical and the spherical geometries to demonstrate the methods. However, the proposed methods work also for problems with other 3D geometries.
The rest of this paper is organized as follows. In the next section we give a brief account of the mathematical models of the drug delivery problems. In Section 3 we propose a discretization method for the diffusion problems based on a finite volume scheme in space and the backward Euler method in time. We then develop a nonlinear least-squares problem in Section 4 for the unknown parameters in the discretized forms to minimize a performance index measuring the 'distance' between the computed drug release profile and a given desirable or experimental release profile. An algorithm will be presented for solving the least-squares problem. We demonstrate the usefulness and accuracy of the numerical methods in Section 5 using various laboratory generated time series of release profiles.
2. Mathematical models of drug delivery. Consider that a device of geometry Ω d loaded with a drug is placed in a container of geometry Ω c filled with a liquid such as water. The drug in the device will diffuse into the liquid. In this section we will first review briefly various mathematical models of drug diffusion from the device into the container. These models, containing unknown parameters, have been detailed in [18,20] and they are summarized below.
i Basic Model (BM): This model consists of only one unknown parameter which is the (constant) effective diffusion coefficient, D. ii Initial Burst Model (IB): In this model, the diffusion process is divided into two phases with two different effective diffusion coefficients: where t c > 0 is the critical time. Besides the two unknown diffusion coefficients D 0 and D 1 , t c is also an unknown parameter in this model. If t c = 0 (or ∞), this model reduces to the basic model. This model is used to handle the socalled Initial Burst Phenomenon in laboratory experiments in which excessive drug left on the surface of a device during the drug loading process causes the diffusion to be much faster in a short initial period of time than in the rest of the process. iii Boundary Layer Model (BL): In this case the liquid region is divided into two subregions. The one around the device, called the boundary layer, is a diffusion dominated region and the subregion outside the boundary layer is convection-dominant so that the drug concentration is uniform in it. If we use Ω l to denote the union of the device and the boundary layer, then the two subregions are respectively Ω l \Ω d and Ω c \Ω l . In this case, we may introduce one unknown parameter to be defined later in this paper. We now consider two different fluid systems: the rotating fluid system and the flow-through system. The mathematical models for these two systems are given as follows.
2.1. Rotating fluid system. A rotating fluid system is a closed fluid system in which the medium is stirred with a magnetic stirrer and not replaced by fresh medium. It is sometimes called a circulatory fluid system. Research on drug release rates in such a system has been conducted in [10,11]. The development of the mathematical models and the estimation of the unknown parameters for this system can be found in [18,20,19,14]. Figure 1 show the device-container settings of cylindrical and spherical devices used in the drug release experiments. The diffusion problem in Cartesian coordinates for a general geometry setting in 3D can be written as where C(x, t) is the drug concentration, V d denotes the volume of the device Ω d , C(t) denote the unknown uniform concentration in the well-mixed region and T is a positive number representing the final time point. This model is similar to the one used in [19] for 2D problems.  2.2. Flow-through system. A flow-though system is one in which the mixed medium is pumped out of the system at a constant rate and fresh medium is pumped into it at the same rate. Therefore, there is a loss of drug mass due to the exchange of mixed and fresh media. Compared to the rotating fluid system, the flow-through system is better related to in vivo [7,17]. Figure 2 depicts a typical flow-through system with a cylindrical device which was used in [16]. The fluid in this system is also stirred and the mass loss is assumed to occur only in the uniform concentration region which is outside the boundary layer. A 2-dimensional drug diffusion problem of this system has been solved numerically in [14]. The drug diffusion in the system is determined by the following initial and boundary value problem where P (t) denotes the rate of drug concentration increase due to the flux across the boundary of Ω l , q(t) denotes the pump rate and the other parameters are the same as those in (2). In this work we assume that q(t) ≥ 0 is either constant or piecewise constant.
3. Discretization. In this section we present a discretization method for the systems given in the previous section. This discretization is based on a finite volume scheme in space and an implicit finite difference scheme in time. We discuss these schemes separately.
3.1. Finite volume discretization. In this subsection we use the finite volume method proposed in [13] for the spatial discretization of the diffusion equations in (2) and (3). In what follows, we will use the cylindrical region depicted in Figures 1  and 2 to demonstrate the finite volume method. We start this discussion with the partitioning of the solution domain.
Since the solution domain Ω c is cylindrical, we partition it into polygonal prisms as follows. Without loss of generality, we assume that the cylinder is parameterized by .., N xy be a set of N xy nodes in x 2 +y 2 ≤ r 2 c for a positive integer N xy . Using this set of nodes, we construct a 2D Delaunay triangulation (or mesh) consisting of triangles and rectangles ( [4]). Then, the interval (0, h c ) in the z-direction is divided into subintervals with N z mesh nodes 0 = z 1 < z 2 < · · · < z Nz = h c . Using the triangulation in the xydirection and the partition in the z-direction we define a partition for the cylinder Ω c so that it contains either triangular or rectangular prisms. This mesh contains N = N xy × N z mesh nodes (x m , y m , z n ) for m = 1, 2, ..., N xy and n = 1, 2, ..., N z , and is referred to as a Delaunay mesh for the solution domain Ω c . We re-order these nodes using a single subscript in such a way that are on the boundary ∂Ω c of Ω c . Note that the union of all the elements of this partition form a computational domain, denoted as Ω h c , approximating Ω c due to the approximation of the curved boundary by facets.
We now define another mesh dual to the Delaunay one defined above. For each mesh node x i := (x m , y m , z n ), we define the Dirichlet tessellation ( [5]) (or Voronoi polyhedron) as follows Clearly, d i consists of points of Ω h c closer to x i than to other mesh nodes. A typical local structure of the Delaunay mesh and Dirichlet tessellations/Voronoi polyhedra is depicted in Figure 3(a) where solid lines represent the Delaunay mesh and the dotted lines represent the corresponding Dirichlet tessellation. Figure 3(b) is an example of the Voronoi polyhedron d i associated with a node x i with 7 neighbouring nodes x j k , k = 1, 2, ..., 7. Each facet of the boundary ∂d i of a tile d i is perpendicular to an edge connecting x i and its neighbour x j . Let e i,j and l i,j denote respectively the edge from x i to x j and the facet of d i perpendicular to e i,j , as illustrated in Figure 3(c).
Applying the one-point quadrature rule to the first term and integrating the second term by parts, we have from the above equality where C i (t) denotes an approximation of C(x i , t), |d i | is the volume of d i and n denotes the unit vector outward-normal to ∂d i . Let I i denotes the index set of all the neighbouring nodes of x i . For each i, ∂d i consists of a finite number of facets l i,j for j ∈ I i , as depicted in Figure 3(b). Let e i,j denote the unit vector from x i to x j . From the construction of the two meshes we have that e i,j coincides with the unit vector outward normal to the facet l i,j . Therefore, the surface integral in (4) can be approximated by the following finite differences where |l i,j | denotes the area of the facet l i,j and |e i,j | denotes the length of edge e i,j . Replacing the surface integral in (4) by the above approximation and dividing both sides of the resulting equation by |d i | , we have for i = 1, 2, ..., N l . This is a semi-discretized form of the diffusion equation.

Time discretization.
We now consider the time discretization of (5). Let Applying the backward Euler scheme on this mesh to (5) gives i defined by the given initial condition. Rearranging the above equation, we have for k = 1, 2, ..., K and i = 1, 2, ..., N l . Equation (6) is a linear system for C k i for i = 1, 2, ..., N l and a fixed k. To solve this system, it is necessary to define boundary conditions. From (2) and (3) we see that C k i =C(t k ) on or outside ∂Ω l (or for i > N l ). Therefore, (6) can be written in the following matrix form: sponding to second term in (6), I is the N l × N l identity matrix and b k is an N l × 1 matrix representing the boundary condition in terms ofC in (2) or (3) which will be discussed in detail below. For the above system, we have the following theorem.
Theorem 3.1. The system matrix of (7) is an M -matrix.
The proof of this theorem can be found in [13]. Because of the above theorem, (7) is uniquely solvable for any k and the discretization satisfies a discrete maximum principle.
3.3. Determination of boundary conditions for (7). In the formulations of (2) and (3), we assumed that the concentration in the region outside the boundary layers ('well-mixed region') is uniform, i.e., C(x, t) =C(t). However,C(t) is unknown and needs to be determined numerically. To determineC(t), we need first to characterize the boundary layer region Ω l \Ω d . For the device-container setting depicted in Figure  4(a), we introduce θ, h l and r l satisfying r d < r l < r c and Therefore, the boundary layer region can be determined by the parameter θ. Using the cylindrical coordinate system centered at the center of the bottom surface of the container (and the device), the domain Ω l in this case is then defined as For the device-container setting depicted in Figure 4(b), we introduce parameters θ and r l satisfying r d < r l < r c and Using this θ and the spherical coordinate system centered at the center of the spherical device, we define the region Ω l as From the formulation of the diffusion problems we have that the concentration outside the boundary layer region Ω l is uniform. Therefore, the drug concentration outside the boundary layer for the rotating fluid system (2) is simply defined as the average concentration over the region given bȳ where |Ω c \Ω l | denotes the volume of Ω c \Ω l as defined before. The concentration C(t) defines the boundary condition for (7). Note that C(x, t) is unknown and thus (8) and (7) are coupled systems for C k andC k :=C(t k ). Recall the assumption that the mesh nodes are ordered in such a way that {x i } N i=N l +1 is the set of nodes outside Ω l . Thus, to determine C k andC k numerically, we decouple (8) and (7) using the following algorithm: withC 0 = 0 defined by the initial condition, and 2. solve The approximation ofC(t) in the region Ω c \ Ω l for the flow-through system can be achieved using a decoupling algorithm similar to Algorithm 3.1. In this case, there is a continuous mass loss due to the pumping at a piecewise constant rate as given in (3d). Applying the forward Euler's finite difference scheme to (3d) , we for k = 1, 2, ..., K − 1 with the initial conditionC 0 = 0. Now we need to define an approximation to P k−1 . Recall that P (t) is the rate of drug concentration increase due to the flux across the the boundary ∂Ω l of the diffusion layer. The exact value of this rate is unknown. However, the term ∆t k P k−1 on the RHS of (11) represents the concentration increase in Ω c \ Ω l from t k−1 to t k , and soC k−1 + △t k P k−1 is the concentration in Ω c \ Ω l at t k without counting the mass loss in the period (t k−1 , t k ). This concentration can be approximated by the RHS of (9) and thus we have approximation from (11) Using Algorithm 3.1 with (9) in Step 1 replaced with (12) denote a given time series of the total mass in the liquid region Ω c \Ω d recorded during a laboratory experiment at the pre-defined experimental time points t e k , k = 1, 2, ..., K e with K e a positive integer. Without loss of generality, we assume that {t e k } Ke k=1 is a subset of the time mesh points {t k } K 1 defined in the previous section. At t e k , using the meshes defined in the previous section we approximate the total mass in the liquid region as follows whereĈ k (u) is the solution from Algorithm 3.1 for a fixed u.
In the rotating fluid system, it is easy to see that the total mass released from the device in infinite time is Therefore, the total mass in (14) can be normalized as follows Similarly, the normalized experimental data are given by In computation, we may approximate M ∞ in (16) by M Ke e when the last laboratory observation time point t e Ke is sufficiently large. Clearly, both R e → 1 and R → 1 as t approaches infinity.
To determine the unknown u, we post the following weighted least-squares problem: for the rotating fluid system Ke k=1 M k e − M k (u) 2 w k for the flow-through system (17) whereĈ k used in R k (respectively in M k ) is determined by (9) and (10) (respectively (12) and (10)), w k 's are weights and u min and u max are given lower and upper bounds on u. Recent advances in numerical methods for solving nonlinear least-squares were discussed in [23]. To solve this least-squares problem, we use the Lavenberg-Marquardt method [9,12]. This method is presented below using the least-squares problem for the rotating fluid system. The cost function in (17) for the rotating fluid system can be written in the following matrix form where Given an approximation u to the exact solution of the least-squares problem, we determine an update δu for u using the following Taylor expansion where || · || denotes the Euclidean norm and J(u) := ∂R(u) ∂u is the K e × N p Jacobian matrix defined by Omitting the term O( δu )δu and substituting the truncated form of (19) into (18), we have The gradient of E with respect to δu is Setting ∇E(u + δu) = 0 and rearranging the resulting equation, we then have which is the Newton's updating formula. Based on this, we define the following Levenberg-Marquardt correction formula [9,12]: where λ ≥ 0 is a constant and I the N p × N p identity matrix. Finally, given an initial guess u 0 to the solution of the least squares problem, for any s = 0, 1, ..., we define the (s + 1)th iterate as u s+1 = u s + αδu s where 0 < α ≤ 1 is a damping parameter and δu s solves (20) with u = u s . The Jacobian matrix J(u) that is needed in the updating formula can be evaluated using divided differences as discussed in [19]. 5. Results and discussion. In this section we use our methods developed in the previous section to estimate effective diffusion coefficients and other unknown parameters using some laboratory generated drug diffusion time series. We divide the numerical experiments into three cases -spherical devices in a rotating fluid system, cylindrical devices in a rotating fluid system and cylindrical devices in a flow-through fluid system. Numerical results will be compared with the results from other methods developed by us. For ease of discussion, we use FULL3D, FULL2D and ANAL to denote, respectively, the full 3D numerical method developed in this paper, the full 2D numerical method developed in [19] and the analytical solutions developed in [20] and [18] for simplified cylindrical and spherical device-container geometries.
In our numerical experiments on FULL3D, we use meshes consisting of triangular and rectangular prisms with the following mesh sizes ∆x = ∆y = r c /2 m , ∆z = h c /2 m , ∆t = 40/2 n (seconds), where m and n are non-negative integers, and r c and h c denote respectively the radius and height in centimeters of the cylindrical container. All the computations for FULL3D were carried out in Fortran 77 double precision under Linux environment.

Case I: Spherical Devices in Cylindrical Containers for the Rotating
Fluid System. For this section, we use the same experimental data as in [18]. This data set contains experimentally observed time series for 6 different devices denoted respectively by S2080-05, S3070-05, S4060-05, S2080-10, S3070-10 and S4060-10, of which the suffixes '-05' and '-10' represents two different drug load levels -5 wt% and 10 wt%. For details of these devices, we refer to [18]. The analytical solutions developed in [18] were based on the assumption that both the device and container are concentric spheres, though the container used in experiments was cylindrical. Using FULL3D, we are able to solve the problem in its real geometries -a spherical device in a cylindrical container.
To make sure that our method and computer code work properly, we have performed a convergence study for Model BM (only one unknown parameter D) using the time series S2080-05. We use δD ≤ 10 −8 as the stopping criterion for the Lavenburg-Marquardt method and the initial guess for D is chosen to be D 0 = 10 −6 . The optimal results for the meshes with (m, n) = (2, 0), (3, 1), (4, 2) and (5,3) in (21) are shown in Table 1. As can be seen, the computed least-squares errors decreases as the mesh sizes decrease. To further study the performance of the numerical method, we consider the convergence of the method in terms of the nonlinear iterations for the mesh with m = 5 and n = 3 in (21). Table 2 shows the convergence history of the update vector δD and the effective diffusion coefficient D and Figure  5 displays the convergence of the pointwise approximation of the release profile for Table 1. Computed diffusion coefficients D and least-squares error E for Data S2080-05 using various meshes.   In order to check the numerical stability of the developed scheme with respect to initial guesses, we solve the least-squares problem for data set S2080-05 using the mesh with (m, n) = (4, 2) in (21) starting from 3 different initial guesses D 0 = 10 −6 , 5 × 10 −6 and 10 −5 . The computed optimal effective diffusion coefficients against the number of iterations are plotted in Figure 6, showing that, starting from the 3 initial guesses, the computed diffusion coefficients all converge to roughly the same point.
We are now ready to present the results from FULL3D and compare them with those from the other methods. All the results computed by FULL3D have been obtained using the mesh with (m, n) = (4, 2) in (21). We first consider the application of FULL3D to Model BM. Table 3 contains the computed diffusion coefficients and the least-squares errors for the 6 experimental times series using FULL3D and ANAL derived in [18]. From the table we see that computed diffusion coefficients and the least-squares errors from both methods are very similar to each other, indicating that both methods give consistent results. Table 4 contains the results using the FULL3D and ANAL for Model IB in which the diffusion coefficient is of the form (1) containing three unknown parameters. In our work, we assume that t c only takes values from the first few experimental observation time points. Again, the results from both methods are similar. From the table we see that the least-squares errors are smaller than those listed in Table  3. Also, the initial burst phenomenon is not obvious for devices with suffix '-05', though there are several noticeable differences in t c for these devices. The initial burst is very prominent for devices with the suffix '-10'. This is because these devices have an initial drug loading twice as much as that of the devices with the suffix '-05'.
Finally, we list the computed results from the two methods for Model BL in Table 5 in which θ is the parameter defined in Section 3 characterizing the width of the boundary layer region. Though there are slight differences in the computed θ from the two methods, the computed D and least-squares errors are very close to each other.

5.2.
Case II: Cylindrical Devices in Cylindrical Containers for the Rotating Fluid System. In [20] and [19] we developed respectively an analytical method (ANAL) and a full numerical method (FULL2D) for the 2D diffusion estimation problem in which the device and container are concentric discs. The case considered in [20,19] is a simplified version of the one depicted in Figure 1(a). In this subsection, we solve the full-scale problem in 3D using FULL3D and compare the results with those from ANAL and FULL2D. For brevity, we will solve the problem using only two sets of experimental data labelled with A3 and A4. The time series A3 and A4 can be found in [19].   Computed diffusion coefficients and the least-squares errors for Model BM using the three methods are listed in Table 6. From the table we see that all the methods have comparable least-squares errors. Unlike the situation for the spherical devices, the estimated diffusion coefficients from FULL3D are smaller than the those from the other two methods. We believe that the results from FULL3D are more accurate because in the 3D formulation (2), as well as in practice, the diffusion occurs in the outward normal directions of both the side and top surfaces of the cylinder, while in the 2D models it occurs only in the outward normal direction of the side surface direction. This shows that the 2D methods over-estimate the effective diffusion coefficients and thus the present 3D method is important for accurately estimating the effective diffusion coefficients for cylindrical devices.
To further show the usefulness of FULL3D, we solve the diffusion coefficient estimation problem in all the three models by FULL3D using the time series A4 and record the results in Table 7. Comparing these results with those in [19, Table  5] we see that the 2-dimensional numerical method FULL2D over-estimates the effective diffusion coefficients.  [15]. The problems have also been solved previously using FULL2D and the results from both methods are listed in Table 8. Similar to Case II, the computed diffusion coefficients from FULL3D are smaller than those from FULL2D because of the reason mentioned in Case II above.
6. Conclusion. We have developed a 3-dimensional numerical technique for estimating the unknown parameters in drug delivery problems. Two fluid systems have been considered: the rotating fluid system and the flow-through system. We have tested the numerical methods using sphere-and cylinder-shaped devices for the rotating fluid system and cylinder-shaped devices for the flow-through system. Results from our 3-dimensional numerical technique have been compared with the those from some existing analytical and 2-dimensional numerical methods. The comparisons show that all the methods give comparable results for spherical drug delivery devices. However, the analytical and 2D numerical methods seem to overestimate the effective diffusion coefficients of cylindrical devices.