Mathematical modeling of diffusion mechanisms for drug delivery can be very useful for accelerating product development and for better understanding the mechanisms that control drug delivery from advanced delivery systems.
The integration of Comsol’s numerical environment with programming in Matlab allows the development of customized and easily integrated solution models. These models easily allow diffusion analysis for drug-releasing systems with desired geometries and constitutive properties.
The following report is an extract from a project carried out for the course of Modeling and Simulation of Physiological Systems, held for Medical Engineering at the University of Rome Tor Vergata.
Table of Contents
Introduction
Drug release systems are becoming more and more popular and in silico simulation can greatly aid in drug selection and manufacturing procedures to provide a specific release profile.
Thanks to significant advances in information technology, mathematical modeling of drug delivery is a field of ever-increasing academic and industrial importance with enormous future potential.
In silico optimization of new drug delivery systems can be expected to significantly increase accuracy and ease of application. Furthermore, computer simulations are likely to become an integral part of future research and development in pharmaceutical technology.
It is therefore essential to have the possibility of carrying out analyzes for the composition, geometry, dimensions and preparation procedure of the various types of administration systems, taking into account the desired route of administration, the dose of the drug and the release profile.
Background
The controlled release of substances originates around 1950 with the development of polymeric systems capable of carrying out a programmed release [1].
Drug-eluting systems
This control can be applied to various fields such as cosmetics, pharmaceutical industry, agriculture and food sciences.
In particular, by focusing on drug delivery systems, these allow for a gradual and continuous release over time, guaranteeing the dose within the optimal ranges. Such systems can be developed in such a way as to optimize the kinetics of the release ensuring the maximum efficacy of the drug therapy.
There are different types of controlled release systems based on different phenomena such as diffusion, dissolution, ion exchange and control with chemical mediators [2].
Therefore, in this report the focus is on diffusion-based systems. Furthermore, such systems are based on the assumption that mass transport limits the drug release rate, the conditions modeled are such for the duration of the release, the device does not over-deflate during release, and the diffusion coefficient remains constant in space and time [3].
Diffusion-convection
Diffusion is the migration of a substance from a region with a high concentration to a region with a lower concentration. The drugs are trapped inside water-insoluble inert capsules consisting of polymeric membranes (reservoir system approach) or within a polymer matrix (monolithic approach). Therefore in both cases the drug diffuses from the charged sphere to the solution.
Therefore the diffusion mechanisms are regulated, in space and time, by Fick’s law. This law links the solute current density to the concentration gradient, by means of the diffusion coefficient:
\begin{equation}
\mathbf{j}_{d}=-D \nabla c
\end{equation}To which is added a convective contribution:
\begin{equation}
\mathbf{j}_{c}=c \mathbf{q}
\end{equation}So by focusing on a control volume we can consider the mass balance. The mass will be the spatial integral of the concentration:
\begin{equation}
m_{p}(t)=\int_{P} c(\mathbf{x}, t) d v
\end{equation}That is, this approach allows you to bring the derivative inside the integral, using the Reynolds theorem.
\begin{equation}
\frac{d m(t)}{d t}=\int_{p} \frac{\partial c(\mathbf{x}, t)}{\partial t} d v=-\int_{\partial P} \mathbf{j} \cdot \hat{n} d a+\int_{P} \mathbf{r} d v
\end{equation}Pertanto la variazione della massa sarà legata ai contributi entranti dal bordo più un eventuale contributo di sorgente. Ovvero, applicando il teorema della divergenza, la conservazione della massa può essere riscritta come:
\begin{equation}
\int_{P}\left(\frac{\partial c}{d t}+\operatorname{div} \mathbf{j}\right) d v=\int_{P} \mathbf{r} d v
\end{equation}Localization
Therefore, by exploiting the arbitrariness of the control volume, this can be localized:
\begin{equation}
\frac{\partial c}{\partial t}+\operatorname{div} \mathbf{j}=\mathbf{r}
\end{equation}And, overall, for an incompressible fluid, a differential problem is obtained:
\begin{equation}
\left\{\begin{array}{l}
\frac{\partial c}{\partial t}-D \Delta_{2} c+q \cdot \nabla \mathbf{c}=\mathbf{r} \\
\text { +B. C. } \\
\text { +I. C. }
\end{array}\right.
\end{equation}The initial conditions, that is the distribution of the species in the system, and the boundary conditions must also be considered. Then this partial differential equation can be easily solved with a finite element solver.
Furthermore, for the description of a possible coating of the volume, a sufficiently thin three-dimensional shell with its own diffusion coefficient could be considered. However, it is more convenient to consider an interface, that is a 2D surface.
A coating admittance applies to the description, for which the following applies:
\begin{equation}
\mathbf{j} \cdot(-\hat{n})=h[[c]]
\end{equation}Often, this value is considered to be:
\begin{equation}
h=\frac{D_{\text {barriera }}}{\text { spessore }}
\end{equation}To which the continuity equation is added:
\begin{equation}
[[\mathbf{j} \cdot(-\hat{n})]]=0
\end{equation}Comsol Multiphysics
This diffusion-convection model is implemented within the COMSOL Multiphysics environment [4].
This environment allows you to perform simulations using advanced numerical methods and using the graphical interface through which to interact by combining different physical phenomena within the same analysis.
It is widely used in science and engineering to simulate and design devices and processes in many application scenarios. In addition to the application of model creation, the integration with the Matlab environment is exploited by transferring the setting and solution of the model to Matlab programming.
Results
To have a model that is easy to set up and possibly also to modify, a block structure is set, each one defined by a subroutine. The block diagram is present in fig. 3.
| Variable | Value | Meaning |
|---|---|---|
| R_cyl | 3E-2 | Cylinder radius |
| H_cyl | 2E-2 | Cylinder height |
| R_sph | 0.3E-2 | Spheres radius |
| D_liq | 1E-6 | Diffusion coefficient for solution |
| D_sph | 2E-6 | Diffusion coefficient for spheres |
| C0 | 3 | Initial condition |
| Y_ctg | 0.1E-3 | Coating admittance |
| d_ctg | R_sph/1000 | Coating thickness |
The code starts by loading the Comsol packages and activating the path to the folder containing the subroutines. After the first initialization phase, the simulation parameters, present in Tab. 1, are set and loaded via mydata.m. These parameters are also reloaded into Comsol via set_parameters(). This, like other precautions below, allows you to re-open this file in Comsol and keep all the values previously saved in Matlab.
Once the parameters have been loaded you can move on to set the geometry. The create_geometry() routine is responsible for creating the geometric model within Comsol. The N spherical geometries ( 𝑁 = 4 in this analysis) and the cylindrical geometry are created considering the parameters (position and dimensions) previously mentioned.
This routine also allows the selection of domains and edges of geometries. Each geometric element created defines an internal domain and some edge elements through a numeric ID which is necessary to apply edge conditions or initial conditions.
Numerical modeling
In particular, the technique to identify them exploits the construction of a bounding box around the single elements and the request, to the Comsol model, of the domains (or edges) inside it.
Finally, this routine returns the updated model with the geometries and IDs of the domains and edges of the spheres and cylinders
Next, the physics is added via the create_physics() routine. The diffusion properties are added to the cylinder and to the spheres, through the diffusion coefficient, considered isotropic, and the initial conditions on the spheres.
So the geometry is meshed with more or less accuracy, chosen through a numeric ID from 1 (Extremerly Fine) to 0 (Extremely Coarse). An ad hoc routine is also implemented to convert the textual flag into the numeric ID.
Time dependencies
The numerical study is then created, the start and end times and the time discretization are set and then the model is solved. The solution routine creates a time dependent study and calls Comsol’s various proprietary numerical solvers.
Finally, the model is passed to the post pocessing routine. The concentration is integrated into the volume to obtain the mass over time, the trend of which is shown in fig. 6a. The concentration in the radial direction in the spheres is also read. To do this, ad hoc cutlines are created that cut the spheres from the center to the edge. Subsequently, the data are extracted exactly on the cutlines and are shown in fig. 6b.
Convergence analysis
A mesh sensitivity analysis is then carried out by repeating, within a loop for, the analysis by varying the degree of refinement of the mesh. The data for each simulation are then collected, such as the number of degrees of freedom, and the final residual mass in the spheres at the last instant of time is considered as a value of interest. The error is then calculated with respect to the assumed true value, that is the simulation data with the densest possible mesh.
From whose results, in fig. 5a and in Tab. 2, the mesh with numerical ID equal to 2 is chosen. An Extra Fine discretization with 154456 degrees of freedom is then chosen (fig. 5b).
| Mesh refinement | Extr. coarse | Extra coarse | Coarser | Coarse | Normal | Fine | Finer | Extra fine | Extr. Fine |
|---|---|---|---|---|---|---|---|---|---|
| ID | 9 | 8 | 7 | 6 | 5 | 4 | 3 | 2 | 1 |
| DoF | 797 | 2520 | 4493 | 8924 | 13044 | 19697 | 45692 | 154456 | 718178 |
| Error | 8.2220 | 1.1591 | 0.5191 | 0.2223 | 0.2215 | 0.2036 | 0.0913 | 0.0214 | 0 |
| Time cost | 5.9820 | 4.5229 | 4.4123 | 5.0982 | 5.8378 | 6.4368 | 10.1138 | 27.1032 | 133.6901 |
This solution has an error of less than 5% compared to the solution considered true but with a computational cost 5 times lower.
Mass trend
Once the numerical simulation has been solved, it is possible to analyze the variation of different quantities.
The trend of the mass within the solute is analyzed, which progressively increases as a function of the drug release by the spheres. The trend is shown in fig. 6a. Then the same quantity is then also read within the four spheres where it progressively decreases in a similar way (fig. 6b).
Therefore, the concentration values are also read, in a radial direction, within the four spheres. The trend is similar to each other and only the trend for sphere number 4 is reported, in fig. 4. The surface color maps are also shown with the concentration at the final and final instant, in fig. 7.
Surface coating
The possibility of inserting a surface coating is then added, i.e. an interface on the sphere that alters its diffusion properties.
This peculiarity, the parameters of which are present in Tab. 1, is directly passed on to routines for the creation of physics indicating the number of spheres (from 1 to 4, with respect to fig. 8a) with respect to which they have been initialized.
Subsequently, the routine enters the coating properties using the thin barrier option and providing Comsol with the related parameters. The user is then returned a visual feedback such as to indicate the spheres on which the coating is applied (fig. 8b).
Coating effects
This coating considerably slows down diffusion and this is evident both from the total mass in the cylinder (fig. 9a) and from the residual mass in the spheres (fig. 9b). That is, the mass in the spheres progressively reduces but does so with different speeds.
Furthermore, analyzing the trend of the concentration in the spheres, in fig. 10 it can be seen how the concentration in the spheres with the coating remains higher over time.
Time extension
Analyzing fig. 9b it is evident how the problem is temporally incomplete so the analysis is extended from 4 to 40 seconds showing how the mass in the two spheres tends to balance itself at the same value, albeit with different speeds. Similarly, the mass in the cylinder tends to a situation of equilibrium (fig. 11).
In appendix are shown the routines for the numerical extraction of the results from the Comsol model.
Conclusions
The integration of Comsol’s numerical environment with programming in Matlab allows the development of customized and easily integrated solution models. These models easily allow diffusion analysis for drug-releasing systems with desired geometries and constitutive properties.
The study of drug release and its diffusion allows to analyze the variation over time and its concentrations in the main solution.
References
- [1] Cong Truc Huynh and Doo Sung Lee. “Controlled Re- lease”. en. In: Encyclopedia of Polymeric Nanoma- terials. Ed. by Shiro Kobayashi and Klaus Müllen. Berlin, Heidelberg: Springer Berlin Heidelberg, 2015, pp. 439–449. ISBN: 978-3-642-29647-5 978-3-642- 29648-2. DOI: 10.1007/978-3-642-29648-2_314. URL: http://link.springer.com/10.1007/978- 3-642-29648-2_314
- [2] Ronald A. Siegel and Michael J. Rathbone. “Overview of Controlled Release Mechanisms”. en. In: Fundamen- tals and Applications of Controlled Release Drug De- livery. Ed. by Juergen Siepmann, Ronald A. Siegel, and Michael J. Rathbone. Boston, MA: Springer US, 2012, pp. 19–43. ISBN: 978-1-4614-0880-2 978-1- 4614-0881-9. DOI: 10.1007/978-1-4614-0881- 9_2. URL: http://link.springer.com/10.1007/ 978-1-4614-0881-9_2
- [3] Juergen Siepmann and Florence Siepmann. “Modeling of diffusion controlled drug delivery”. en. In: Journal of Controlled Release 161.2 (July 2012), pp. 351–362. ISSN: 01683659. DOI: 10 . 1016 / j . jconrel . 2011 . 10.006. URL: https://linkinghub.elsevier. com/retrieve/pii/S0168365911009588 .
- [4] Comsol. “Combining Convection and Diffusion Ef- fects”. In: (2022). URL: https : / / www . comsol . com / multiphysics / convection – diffusion – equation.
Appendix
Coating
The coating is applied only if the variable 𝚌𝚘𝚊𝚝𝚒𝚗𝚐 contains the indices of the spheres on which to apply the Thin Diffusion Barrier interface.
if not(sum(coating==0)) surface_with_coat =[];
for i=1:length(coating)
surface_with_coat=[surface_with_coat , sph_boundary_alt(:,coating(i))’];
end
model.component(’comp1’).physics(’tds’).create(’tdb1’, ’ThinDiffusionBarrier’, 2);
model.component(’comp1’).physics(’tds’).feature(’tdb1’).set(’ds’, num2str(d_ctg));
model.component(’comp1’).physics(’tds’).feature(’tdb1’).setIndex(’Ds’, D_ctg, 0);
model.component(’comp1’).physics(’tds’).feature(’tdb1’).selection.set(surface_with_coat); % graphics to hihglight coated spheres
figure;
mphviewselection(model,’geom1’,surface_with_coat ,’boundary’,’facealpha’,0.5,’
facecolorselected’,[0 1 0]) title(’Coated sphere’)
endIntegral of concentration
The extraction of the mass through the integration of the concentration requires the mphint2 routine.
mass =mphint2(model,’c’,’volume’,’selection’,1);
figure; plot(initValue:step:finalValue ,mass);
title(’Mass over time’); xlabel(’Time [s]’); ylabel(’Mass [mol]’)Cutline
The extraction of the concentration along a line requires the construction of the line and the extrapolation of the results on it via mphinterp.
for i=1:length(C_sph(1,:))
% create cutline to read the data
cutLine=strcat(’cln’,num2str(i));
model.result.dataset.create(cutLine, ’CutLine3D’);
% syntax ’genpoints ’, value , point index , coordinate
model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(1,i), 0, 0); model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(2,i), 0, 1);
model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(3,i), 0, 2); model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(1,i)+R_sph, 1, 0);
model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(2,i), 1, 1); model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(3,i), 1, 2);
% extract results
[conc, space,time_plot,unit]=mphinterp(model,{’c’,strcat(’cln’,num2str(i),’x’),’t’},’dataset’
,cutLine);
[space,index]=sort(space’,’ascend’); % the data are rearranged to have the correct x distance
x_to_plot=space; y_to_plot=conc’; % and the same index is used to rearrange the corresponding concentration values
for j=1:length(y_to_plot(1,:))
y_to_plot(:,j)=y_to_plot(index(:,j),j);
end
figure; plot(x_to_plot ,y_to_plot)
xlabel(unit{2}); ylabel(unit{1}); legend(string(time_plot(:,1)’),’Location’,’eastoutside’) title(strcat(’Radial concentration for sphere n.’,num2str(i)))
end
