Numerical prediction of drug release from polyurethane samples combining kinetic model developed by experimental data

With the aim of optimizing Drug Eluting Stents (DES), particular attention has been laid on computational methods of controlling the drug release profile among researchers. Consequently, various models and simulations are available in the literature. Nevertheless, validations based on biorelevant in-vitro trials are lacking. In the present study, a comparison of drug release from polyurethane samples between calculated results and experimental-data has been carried out. The calculation results are from a numerical simulation and a newly established mathematical model for reproducing the liberation kinetic. Different fluid flow rates and initial drug concentrations in polymer have been taken into account. ⁎ Corresponding author: jianfei.song@ensam.eu


Introduction
Controlling the drug release profile from DES acts a crucial role to improve DES performances [1]. The numerical and experimental works have been widely conducted focusing on the drug release behaviors from DES in literature [2][3][4]. However, the validations by biorelevant in-vitro tests are still scarce because of the limited applicable mathematical models to depict the drug delivery manner. Even though a number of empirical models have been put forward through experimental data for a range of drug carriers, the improvement of these models is necessarily required as they all remain certain limits to a more widespread application [5][6][7]. Within the limited validation works in literature, P. R. S. Vijayaratnam et al. [8] compared the drug distribution inside drug carrier and tissue between simulation and experiment in 2018. Nevertheless, the discrepancy of built numerical model from the experimental rig makes the conclusion less convincing. In addition, in 2019 C. M. McKittrick et al. [9] achieved prediction of the specific diffusion coefficient within drug carrier via numerically fitting the tested drug release data. Accordingly, a good validation has been achieved for drug concentration inside tissue between in-vivo and numerical results by means of the fitted diffusion coefficient. However, a predictable kinetic model of drug release is still absent to precisely depict the drug release profile without relying on the data set. So as to achieve this, the fundamental mechanisms and related influencing factors should be better understood which makes it further difficult. Within the recently published paper by N. Abbasnezhad et al. in 2020 [10], a new predictive kinetic model has been proposed to describe the drug release profile. Therefore, based on this mathematical model for reproducing drug release process, a comparison between calculated results and test-data of a model of drug released from polyurethane samples has been carried out in the present study. The calculation data is from this drug liberation kinetic model along with the numerical simulation. Concerning the influential factors, different fluid flow rates and initial drug concentrations have been taken into account.

Test rig and material
The experimental facility adopted is displayed in Figure 1 with several parallel test tubes. The test tube is fabricated in rectangular shape with length of 130 mm and the squared cross section with length of 30 mm.

Fig. 1. Test rig.
Fluid property: Water was adopted as the fluid flow in these tests to take away the medications from polymer. The density and dynamic viscosity are 1000 kg/m 3 and 6.9e-4 Pa.s respectively. The flow rate has been set to be 0 ml/s, 6.5 ml/s, 7.5 ml/s and 23.5 ml/s separately.
Polymer material: Non-degradable polyurethane (PU) was employed as the drugs carrier. Polymer samples with dimensions of 30×5×2 mm3 have been used with a certain drug dosage defined as mass ratio of drug/(drug+polymer): 10% and 20% respectively. The drugs have been homogeneously distributed in polymer fixed at the center of flow channel whilst contacting the bottom surface.
Choice of drug: The drugs loaded in polymer were selected to be diclofenac epolamine (C20H24Cl2N2O3). In addition, the drug particle is in granular shape with density of 450.7 mg/ml. The drug solubility at temperature of 37 o C in water is around 5.554 g/L.
Proposed kinetic model: The proposed kinetic model in [10] is described in equation (1): Where t M denotes the drug release mass from polymer which is calculated by the mass difference between initial and left drug dosage in polymer along with tile evolution. 0 M represents the initial drug mass in polymer.  is defined as

95
T is defined as follows: indicates the initial drug amount in polymer.
Re describes a reduced Reynolds: More details regarding this kinetic model can refer to literature [10]. Additionally, the drug release governed by diffusion can be expressed as: Where Sec is the hydraulic surfaces of polymer. Therefore, the diffusion coefficient D can be depicted as: Sec is well defined to be  (6): Accordingly, the predicted diffusion coefficients under different flow rates and different initial drug dosages are listed in Table 1 complying the available experimental data:    Figure 3 shows the released drug mass percentage from polymer at t=0.1 d versus these different mesh sizes when the flow rate is set to be 7.5 ml/s and the initial drug dosage is 10% in polymer. As observed, when the mesh elements is over 8.6 million, the accuracy can be achieved within 1%. Therefore, the mesh size of 8.6 million is finally determined with mesh distribution displayed in Figure 4. To be coherent with experiments, water with dynamic viscosity of 6.9e-4 Pa.s and density of 1000 kg/m 3 has been considered in numerical flow domain as it is regarded to be incompressible and Newtonian fluid. Navier-Stokes equations including mass conservation and momentum equations are adopted to govern the fluid flow as described in equation (7) and (8). Different flow rates of 0 ml/s, 6.5 ml/s, 7.5 ml/s and 23.5 ml/s have been all tested corresponding to the inlet sectional mean velocity of 0 m/s, 0.0072 m/s, 0.0083 m/s and 0.026 m/s respectively. The relevant Re number is below 1200 within laminar regime. An uniform velocity profile is considered at inlet and the polymer is fixed at a region where the flow has been fully developed. A constant pressure is set at the outlet and the artery wall is deemed to be rigid.
Where V indicates the flow velocity, p denotes the pressure,  represents the flow density and  is the dynamic viscosity. Specific to the drug transport in lumen, diffusionconvection equation has been employed as shown in equation (9) [11]:  As for the boundary conditions, the drug concentration is fixed to be 0 at the inlet, while a drug flux of 0 is set at the outlet. At the interfaces between polymer and flow domain, the drug flux is thought to be continuous. Initial drug concentrations in polymer are set to be 202 mol/m 3 and 356 mol/m 3 corresponding to 10% and 20% separately. In the flow domain, the initial drug concentration is set to be 0. Comsol 5.1 with finite element method has been adopted for modeling, meshing and computing. The iterative methods have selected generalized minimal residual method (GMRES). The streamline upwinding Petrov-Galerkin (SUPG) scheme and the backward differential formula (BDF) scheme with variable order are imposed for achieving the spatial and temporal discretization respectively of the governing equations. The residual value is set to be 10 -5 for both continuity and momentum equations. The residual value of 10 -4 has been imposed for drug diffusion. The initial time step is firstly set to be 0.001 s which value can be adaptively adjusted during computing process counting on the physics and scheme employed.

Results and discussion
The flow topology around the polymer has been investigated in Figure 5 under different flow rates. A slice located at the middle of flow channel has been made in the flow direction. As observed that the vortices have been developed both upstream and downstream the polymer as the flow has been strongly disturbed by this obstacle. In addition, the increased flow rate is prone to bring larger recirculation region.   Similarly, Figure 7 shows the comparison of drug release profile from polymer among the numerical, experimental and kinetic data with different flow rates and initial drug dosage of 20%. Still, a good validation of drug release from polymer is achieved between the experimental and calculated data at flow rates of 0 ml/s, 6.5 ml/s, 7.5 ml/s and 23.5 ml/s respectively. Increased flow rate boosts the drug release process with decreased drug release time. Specific to the drug release period with the same flow rate compared to that in Figure 6, it is observed that the augmented initial drug dosage brings decreased drug release period. The reason is that higher drug dosage can increase the pores inside polymer structure as proved in [3].

Conclusions
In the present study, drug release profiles from polymer have been focused and compared between calculated and experimental results considering different influential factors: different flow rates and different initial drug dosages. A good validation of drug release from polymer has been well achieved at all the conditions concerned based on the predictive kinetic model of drug release reproduction put forward by us previously. It is also proved that the increased flow rate and initial drug dosage can reduce the drug release period. Moreover, disturbed flow has been observed both upstream and downstream the polymer. The recirculation region tends to be increased by the larger flow rate.