Mathematical Model for Estimating Parameters of Swelling Drug Delivery Devices in a Two-Phase Release

Various swelling drug delivery devices are promising materials for control drug delivery because of their ability to swell and release entrapped therapeutics, in response to physiological stimuli. Previously, many mathematical models have been developed to predict the mechanism of drug release from a swelling device. However, some of these models do not consider the changes in diffusion behaviour as the device swells. Therefore, we used a two-phase approach to simplify the mathematical model considering the effect of swelling on the diffusion coefficient. We began by defining a moving boundary problem to consider the swelling process. Landau transformation was used for mitigating the moving boundary problem. The transformed problem was analytically solved using the separation of variables method. Further, the analytical solution was extended to include the drug release in two phases where each phase has distinct diffusion coefficient and continuity condition was applied. The newly developed model was validated by the experimental data of bacterial cellulose hydrogels using the LSQCURVEFIT function in MATLAB. The numerical test showed that the new model exhibited notable improvement in curve fitting, and it was observed that the initial effective diffusion coefficient of the swelling device was lower than the later effective diffusion coefficient.


Introduction
The safe and effective delivery of therapeutics at the target site of action is a challenging task for drug delivery scientists. High dose or dosage frequency, undesirable variations in drug concentration, nonspecific distribution of drugs, and severe side effects are the major limitations associated with conventional drug administration. These limitations lead to poor patient compliance and therapeutic response [1][2][3]. To overcome these limitations, controlled release drug delivery devices, such as those using hydrogels, nanoparticles, micelles, liposomes, and membranes, have gained considerable attention in the last four decades [4]. These controlled release systems are designed to deliver drugs at the desired rate to the target site; thus, they can reduce the dosage and toxicity of a drug and enhance drug efficacy [2]. Stimuli-responsive materials that exhibit conformational changes in their structure in response to physical or physiological stimuli, such as pH, temperature, ionic strength, biomolecules The research work on mathematical modelling of controlled drug delivery has been significantly advanced; however, the effect of swelling on the diffusion coefficient is usually neglected in the mathematical model. The diffusion coefficient is a significant parameter in estimating the drug release, especially in a diffusion-controlled release device as the diffusion rate depends on this parameter [23]. The diffusion coefficient is also known as the diffusivity or the ability of a substance to diffuse. The diffusion coefficient can be in a form of a constant or a function of time. Most researchers tend to consider the diffusion coefficient as a constant parameter [17][18][19][20][21]. However, for a device with a significantly high swelling rate, a constant diffusion coefficient might not relevant. A study by Fu et al. [24] showed that, as the swelling ratio of a hydrogel increases, the diffusion ratio of the sample also increases. Theoretically, as the device swells, the pore size gets larger, facilitating the drug diffusion and, eventually, causing an increment in the diffusion coefficient. The study by Ahmad et al. [9] showed that the pore size of a higher swelling device is bigger than devices with lower swelling. The diffusion-porosity relation computed by Ray et al. [25] shows that the diffusion differs with different porosity values and geometries.
Therefore, it is worth developing a mathematical model that acknowledges the dynamic behaviour of drug diffusion as the drug device swells. However, in some conditions, considering the diffusion coefficient as a function of time will not give an analytical solution [26]. Previously, to study the dynamic behaviour of drug diffusion during the initial burst release, Wang et al. [27] divided the whole drug release process into two phases with distinct constant diffusion coefficient in each phase. The first phase is the initial burst release phase and the second phase is the normal release phase. Therefore, in this study, we developed a mathematical model for estimating drug release from a swelling hydrogel-based drug delivery device considering the dynamical constant diffusion coefficient by adopting the idea of two-phase drug release. We solved the fractional drug release model using the advection-diffusion equation, an approach similar to that used by Bierbrauer [20]. The solved fractional drug release model by Bierbrauer [20] was introduced as a one-phase model in this study. Further, the model was extended to include the two-phase release by splitting the drug release process into two phases, as suggested by Wang et al. [27]. Both phases have constant diffusion coefficients with different values.

Mathematical Solution for Drug Release from a Swelling Device
The drug release from a swelling device to an external fluid is influenced by many factors, such as drug diffusion and device swelling processes. The diffusion and swelling of a cylindrical drug delivery device (hydrogel disc) in a 3D space are illustrated in Figure 1.
In this section, two fractional drug release models for swelling devices are developed: one-phase and two-phase models. The one-phase model represents a drug release with a single constant diffusion coefficient throughout the entire release process while the two-phase model represents a drug release with a piecewise diffusion coefficient.

One-Phase Model
We start from a basic model considering a single constant diffusion coefficient. The advection-diffusion equation is a well-known model when considering the drug diffusion and device swelling problem simultaneously. The general advection-diffusion equation is expressed as where ∂c ∂t is the local rate change of concentration over time, c(x, y, z, t) is the concentration of drug in the device, D(t) is diffusion coefficient, c is the concentration gradient, and u(x, y, z, t) is the device swelling velocity. To solve this equation analytically, the diffusion coefficient is assumed as a constant, where D ≡ 0. Thus, For simplicity, the one-dimensional (1D) case with a sink boundary condition is considered. The 1D region of device swelling and drug diffusion is shown in Figure 2. The 1D case model of drug release from a swelling device is then expressed as with a constant initial condition c(x, 0) = c 0 in 0 < x < X(t). No flux boundary conditions were assumed at the left edge (hydrogel centre), ∂c ∂x (0, t) = 0, and the sink boundary condition was assumed at the right edge (hydrogel surface/boundary), c(X(t), t) = 0. The boundary X(t) is moving only in positive x direction starting from X(0) = L. The device swelling is assumed to be uniform.
As the device swells, the velocity gradient is located across the region of the device. The velocity gradient is the difference in the velocities of growth (growing) and initial boundaries. This relation can be expressed as Therefore, the velocity, u is expressed as Using Equation (3), Equation (2) can be rewritten as with the following initial and boundary conditions: As the moving boundary domain of the above model is difficult to solve, Landau transformation is applied in the model, wherein the new parameters are defined as , Introducing the parameter ζ not only helps in handling the original moving boundary problem but also eliminates the advection term, thereby providing a simpler model to represent the drug transport in a swelling device.
with the following initial and boundary conditions: After the transformation, the spatial domain of this model is a constant between 0 and 1, instead of a moving domain that depends on time. Equation (5) can be solved using the separation of variables method. The solution is assumed to be a product of the functions of ζ and τ, From the separation of variables method, two sets of ordinary differential equations (ODE) are obtained: Equations (6) and (7) are solved subject to the boundary conditions, giving Further, the solutions of A(ζ) and B(τ) are rearranged using superposition to get an equation of c(ζ, τ).
The initial condition, c(ζ, 0) = c 0 , is then employed in the above equation, thereby resulting in Substituting the solution of d n into Equation (8), the following equation is obtained: The derived solution, Equation (9), is the equation of the concentration of drugs in the hydrogel. Further, an equation for fractional drug release is derived to study the release rate in a swelling device. Fractional drug release can be defined by where M(τ) is the total mass released and M(0) is the total initial drug loaded. Total mass released, M(τ), can be calculated by subtracting the mass of drug remained in the hydrogel (at any time τ) from the total initial drug loaded, In original variables, this equation can be defined as The integral term in Equation (11) includes the growth boundary function, X(t), which represents the swelling behaviour. In this study, we consider the logistic growth function where m = lim t→∞ X(t) L , L is the initial length of the drug device, and g is the growth parameter. This equation is then integrated and substituted into the fractional drug release equation, expressed in Equation (11). By solving the integration of the logistic function, the fractional drug release, Equation (11), can be expressed as where The solution for the fractional drug release equation, Equation (13), represents the one-phase model.

Two-Phase Model
Now, we extend the one-phase model to the two-phase model. The two-phase model is important to acknowledge the fact that the diffusion coefficient for a swelling drug device cannot remain constant throughout the process. The changes in the device structure due to swelling will affect the drug diffusivity through the device. Theoretically, as the device swells, the device porosity increases, and the diffusion coefficient is also increased. Instead of considering the general D = D(t), our approach is to divide the whole release process into two phases. The constant D differs in each phase where D 0 and D 1 are the diffusion coefficients in the first and second phases. The critical time, t c , is the borderline between the first and second phases of the entire release profiles. In developing this two-phase drug release model, the growth parameter, g, is assumed to be constant throughout the entire release process. As similar methods from solving the one-phase model are to be adapted in this two-phase model, we begin by applying the transformed advection-diffusion equation, Equation (5), in two different phases with two different parameters of diffusion coefficients.
with the following initial and boundary conditions: This model is subject to the initial and boundary conditions, as expressed in Equation (5), with the addition of a continuity condition, c(ζ, τ − c ) = c(ζ, τ + c ). The continuity condition is introduced to ensure a continuous drug release from both phases at t c . We assumed that the solution for the first phase of the model is given by the solution of the one-phase model, Equation (13) with D = D 0 at 0 < t < t c . Thus, we proceed to solve the second phase model, Equation (14b).
Again, the separation of variables method is used to solve Equation (14b). Two sets of solution were obtained for the ODE by applying the fundamental steps for the separation of variables method, which are The solution for the space-dependent function, A(ζ), is the same as that discussed in Section 2.1 because the conditions are the same throughout the device, 0 < ζ < 1. However, the solution for the time-dependent function, B(τ), differs because we consider drug diffusion at the second phase, which begins at τ = τ c instead of at τ = 0. Then, by using superposition, we obtain the following solution: To determine h n , we apply the continuity condition at τ = τ c , and the coefficients obtained from both equations are matched. From the continuity condition, c(ζ, τ − c ) and c(ζ, τ + c ) represent the drugs' concentration in the hydrogel at τ c in the first and second phases, respectively. Function c(ζ, τ − c ) is taken from Equation (9) with τ = τ c and D = D 0 while function c(ζ, τ + c ) is taken from Equation (15) The coefficient h n obtained is The coefficient h n is then substituted into Equation (15), giving This solution gives an equation of concentration of drugs in the hydrogel in the second phase. Further, we substitute Equation (16) into Equation (10) to solve for the fractional drug release in a swelling device after the critical time.
In the original variable, the function for fractional drug release on the second phase is Similarly, for this two-phase model, we consider the logistic growth function, Equation (12), to represent the swelling function, X(t). Therefore, the integral terms in Equation (17) can be further specified as To simplify the equations, let where Thus, Equation (17) can then be expressed as Finally, placing the solutions of fractional drug release from both phases together gives Therefore, the newly established two-phase model, which considers the dynamic diffusion coefficient is complete.

Algorithm for the Two-Phase Model
To determine the growth parameter, g, diffusion coefficients, D 0 and D 1 , and critical time, t c , we implement the least-squares method in the two-phase model. The numerical algorithm for the method is as follows: Step 1 Input the experimental data with discrete value i = 1, 2, 3, ....., K.
Step 2 Find the determined growth parameter, g * , such that it minimises where functionX(t j ) represents the devices' radius at time t j and X(t j , g) is the logistic function defined in Equation (12) using growth parameter g at time t j .
Step 3 Output the determined growth parameter, g * . Use g = g * for the next step.
Step 4 Find (D * 0 , D * 1 ) such that it minimises whereR(t k ) represents the fractional drug release data and R N is where Step 5 Output the optimal parameters, ε opt , D opt , and t c , and stop.

Drug Loading
Bovine serum albumin (BSA), a readily water-soluble model protein drug, was physically loaded (entrapped) into the hydrogel discs by the swelling diffusion method, as described earlier [9]. Briefly, discs were soaked into 25 mL of BSA solution for 24 h at 37 • C. The hydrogel discs were then rinsed with 0.1 M HCl and dried at 25 • C. The concentration of BSA in the loading solution before and after the experiment was determined using Bradford reagent assay. The entrapment efficiency (%) of the BSA in hydrogels was calculated to be 54.5%, 44.4%, and 39.7% for 208035, 307035, and 406035 hydrogels, respectively.

Drug Release
In vitro release experiments were performed by placing the hydrogels into the simulated gastric fluid (SGF) with pH 1.2 for the first 2 h, and then immediately transferred into the simulated intestinal fluid (SIF) with pH 6.8 for the remaining time until maximum drug release. This experiment was done at 37 • C with constant shaking at 50 rpm. The growth of hydrogel and the concentration of the drug in both fluids were measured on an hourly basis. After each hour, 20 mL samples were drawn from release media and Bradford assay was performed to determine the amount of BSA released.

Results and Discussion
The numerical test of the newly developed mathematical model was performed on the experimental data of the three different formulations of hydrogels. The results of our previous study [9] suggest that the drug release in the SGF was negligible because of the lower swelling of hydrogels in acidic pH.
Hence, only the experimental data in SIF were employed for numerical analysis of the developed models. The diffusion and swelling processes are assumed to be uniform and occur only in the radial direction. Thus, the 1D region that we consider is from the centre of the disc (x = 0) to the tip of the device's radius (x = X(t)), as shown in Figure 2.
Both of the newly developed fractional drug release models, that is the one-phase model, Equation (13), and the two-phase model, Equation (19), were tested with the experimental data of the three different formulations of swelling hydrogel-based delivery devices. Parameters g, D 0 , D 1 , and t c were estimated for each hydrogel using the LSQCURVEFIT function in MATLAB software.

Model Test for Device Swelling
First, the growth parameter, g, for all hydrogels was determined using the logistic growth function. The logistic growth function, X(t), represents the radius of the drug device at time t. The logistic function shows an increase in the radius of the drug devices from the original length until it reaches the maximum length. Figure 3 shows that the swelling behaviour of all devices is well-fitted with the logistic function. As presented in Table 1, the determined growth parameter value of the 208035 hydrogel is the highest of all tested devices. This result shows that the 208035 hydrogel swells faster than the other two devices. Meanwhile, the determined growth parameters of the other two hydrogels (307035 and 406035) were similar.

Model Test on Drug Release Profiles
In this section, the one-phase model and the two-phase models are fitted to the experimental data, and the results are compared. Parameters D 0 , D 1 , and t c were determined in this process. Figure 4 shows that, for all tested hydrogels, the curve-fitting of the extended two-phase model fits much better than the basic one-phase model. The diffusion coefficients of both models are determined and presented in Table 2. These results suggest that the determined diffusion coefficient in both models (D, D 0 , and D 1 ) have a descending order for the 208035, 307035, and 406035 hydrogels, thereby indicating that the drug release is the fastest in the 208035 hydrogel among the three tested hydrogels. This result is closely related to the pore size of the tested hydrogels. According to our previous findings [9], the pore size of the 208035 hydrogel is highest (i.e., hydrogels are less dense) as compared to the other two hydrogels where the pore size of 208035, 307035, and 406035 hydrogels are around 50-100, 12-27, and 4-14 µm, respectively.  Table 2 presents that the determined diffusion coefficients of all tested hydrogels are lower in the first phase than in the second phase (D 0 < D 1 ), where D 0 is about 70-90% slower than D 1 . This condition is because of the increase in the pore size of the devices as they swell, thereby resulting in easier and faster diffusion. Furthermore, the determined diffusion coefficient D of the one-phase model is in between the values of diffusion coefficients D 0 and D 1 of the two-phase model. Although the 208035 hydrogel has lower diffusivity in the first phase (D 0 < D 1 ), it appears to release drugs more quickly in the initial phase. As shown in Figure 4b, the steeper slope in the first phase of the two-phase model curve indicates a higher initial fractional drug release rate from the device. The fractional drug release rate can be calculated by where R is the fractional drug release, t is the time, and h is the step size. If the fractional drug release rate in the first phase is higher than the second phase ( ∆R ∆t | t<t c > ∆R ∆t | t≥t c ), we can say that the drug device releases drugs more quickly initially. As presented in Table 3, the measured value of ∆R ∆t of the 208035 hydrogel is higher in the first phase than the second phase ( ∆R ∆t | t<t c > ∆R ∆t | t≥t c ), thereby demonstrating the fast initial drug release. Since D 0 < D 1 for 208035 hydrogel, it shows that this condition is not caused by the diffusivity.
Based on the results in Table 1, we found that the 208035 hydrogel gives the highest value of the determined growth parameter. In our previous study [9], we observed that hydrogels with greater swelling exhibited higher entrapment efficiency, which means that hydrogel with a higher growth rate has more drugs loaded at the beginning. The high drug load might lead to a fast initial drug release rate.
Nevertheless, the resulted curve fitting of the two-phase model improved significantly as compared to the basic one-phase model. This finding is supported by the error values listed in Table 4. The 208035 hydrogel showed the least error improvement as compared to the other two devices. The least-squares error of the two-phase model for the 307035 and 406035 hydrogels were almost ten times smaller than that of the one-phase model. Nonetheless, this new two-phase model manages to reduce the least-squares error for all the tested drug devices.

Conclusions
In this study, we successfully developed a new two-phase drug release model considering device swelling. The performance of the developed two-phase model was numerically tested using the experimental data of three different BC-g-P(AA) hydrogels. It is shown that the newly developed two-phase model exhibited significant improvement to the previously known model of drug release in a swelling device. The least-squares error was reduced, and the established model proved to be capable of estimating the drug release in a swelling drug delivery device. The numerical results also detected that the initial diffusion coefficients for all three hydrogels were lower than the later diffusion coefficients. The smaller diffusion coefficient at the initial phase as compared to the later phase was mainly related to the changes in the device pore size because of swelling. A detected initial burst for device 208035 is believed to be caused by the high growth parameter of that device, which leads to a high amount of drug loaded. The parameters of this model helped in analysing and predicting the behaviour of a swelling drug delivery device. Thus, this model will facilitate the researcher in designing controlled drug delivery devices in the future. This model can be used to determine the parameters of drug release kinetics from a swellable matrix system using the in vitro results.