Next Article in Journal
Lie Symmetry Group, Invariant Subspace, and Conservation Law for the Time-Fractional Derivative Nonlinear Schrödinger Equation
Next Article in Special Issue
Traveling Band Solutions in a System Modeling Hunting Cooperation
Previous Article in Journal
Thermodynamic Modelling of Transcriptional Control: A Sensitivity Analysis
Previous Article in Special Issue
Transport Phenomena in Excitable Systems: Existence of Bounded Solutions and Absorbing Sets
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Modeling of Recursive Drug Delivery with Diffusion, Equilibrium, and Convection Coupling

by
Rosaura Hernandez-Montelongo
1,
Javiera Salazar-Araya
2,†,
Jacobo Hernandez-Montelongo
2 and
Juan Paulo Garcia-Sandoval
3,*
1
Department of Translational Bioengineering, University of Guadalajara, Guadalajara 44430, Mexico
2
Department of Physical and Mathematical Sciences, Catholic University of Temuco, Temuco 4813302, Chile
3
Department of Chemical Engineering, University of Guadalajara, Guadalajara 44430, Mexico
*
Author to whom correspondence should be addressed.
Current address: Higher Education Access Program, University of La Frontera, Temuco 4811230, Chile.
Submission received: 24 May 2022 / Revised: 16 June 2022 / Accepted: 18 June 2022 / Published: 22 June 2022
(This article belongs to the Special Issue Transport Phenomena Equations: Modelling and Applications)

Abstract

:
In this work, a mathematical model to describe drug delivery from polymer coatings on implants is proposed. Release predictability is useful for development and understanding of drug release mechanisms from controlled delivery systems. The proposed model considers a unidirectional recursive diffusion process which follows Fick’s second law while considering the convective phenomena from the polymer matrix to the liquid where the drug is delivered and the polymer–liquid drug distribution equilibrium. The resulting model is solved using Laplace transformation for two scenarios: (1) a constant initial condition for a single drug delivery experiment; and (2) a recursive delivery process where the liquid medium is replaced with fresh liquid after a fixed period of time, causing a stepped delivery rate. Finally, the proposed model is validated with experimental data.

1. Introduction

Polymer coatings on implants and other biomedical devices are frequently used as local drug delivery systems, which allows for in situ release in the immediate targeted tissue, optimizing the therapeutic effects. Medications can be released at controlled rates in the range of effective treatment levels, increasing their bioavailability [1]. In this sense, mathematical modeling of drug delivery and release predictability is useful in development and understanding of drug release mechanisms from controlled delivery systems. In most cases, drug diffusional mass transport is the predominant step in the control of drug release, while in others it is present in combination with polymer swelling or polymer erosion [2]. In general, it is common to use power law equations for modeling release kinetics, such as the Higuchi and Ritger–Peppas models.
The Higuchi model [3] implies various assumptions, namely, that the initial drug concentration in the polymer is higher than the drug’s solubility, that drug diffusion occurs in only one dimension (i.e., the edge effect must be negligible), that the drug particles are smaller than the polymer thickness, that matrix swelling and dissolution are negligible, that drug diffusivity is constant, and that perfect sink conditions are always attained. The Higuchi equation is provided by
M t / M = k t 1 / 2
where the ratio M t / M is the fraction of drug released at time t and k is the Higuchi dissolution constant, which represents a Fickian diffusion of drugs without the matrix dissolution.
The Ritger–Peppas model [4] is a simple relationship to describe drug release from a polymeric system. This model analyzes both Fickian and non-Fickian release of drug from swelling as well as non-swelling materials; however, it is only applicable up to 60% of the drug amount released. The Ritger–Peppas release equation is provided by
M t / M = k t n
where the ratio M t / M is the fraction of drug released at time t, k is the Ritger–Peppas kinetic constant, which characterizes the drug–matrix system, and n is the exponent that indicates the drug release mechanism.
These semi-empirical models are easy to use and the established empirical rules help to explain transport mechanisms. However, they do not provide further insights into a more complex transport mechanism. Furthermore, these models might fail when specific experimental methodology or physicochemical processes are involved [5]. For this reason, it is important to develop a mathematical model to adequately study the release kinetics from polymer coatings to liquid media. In this regard, this work is about a mathematical model that describes drug delivery considering a unidirectional iterative diffusion process in Fick’s second law while considering the convective phenomena from the polymer matrix to the liquid where the drug is delivered and the polymer–liquid drug distribution equilibrium. The resulting model is solved using Laplace transformation for two scenarios: (1) a constant initial condition for a single drug delivery experiment; and (2) a recursive delivery process where the liquid medium is replaced with fresh liquid after a fixed period of time, causing a stepped delivery rate. Finally, the proposed model is validated with experimental data.

2. Model Formulation

The main goal of the following model is to understand the drug release mechanisms from controlled delivery systems. Thus, the model must be able to evaluate whether the delivery system is controlled by the diffusion of the drug in the solid matrix (usually a polymer), by convection from the polymer surface to the liquid, or by the solid–liquid drug distribution equilibrium. In particular, we are interested in the experimental configuration shown in Figure 1, where a polymer coating (solid matrix in the following) on an implant initially impregnated with the drug is immersed in a volume of liquid, V l , for a fixed period of time, producing a time-varied drug concentration in the liquid bulk, C l t , that can be measured.
In order to model drug delivery, the following assumptions are considered:
  • Drug diffusion is only in one dimension, x (i.e., the edge effect must be negligible); therefore, the concentration of the drug in the solid matrix is expressed as C x , t
  • Drug particles (solute) are smaller than solid thickness
  • Matrix swelling and dissolution are negligible, therefore, the length of the diffusion region, 2 L , is constant, and by symmetry the problem is analyzed in the following domain: 0 x L
  • Drug diffusivity in the solid matrix, D, is constant and Fick’s law can be applied, i.e., J = D C / x
  • There exists an equilibrium between the concentration of the drug in the interface of the solid matrix and the liquid, provided by C e q = γ C l
  • There is a convective rate that depends on the difference of the concentration in the solid–liquid interface, C x = L , and the equilibrium concentration, C e q , with the form J x = L = k C C x = L C e q
  • The entire surface of the solid is completely submerged.
In the following sections, a model is developed for two scenarios: (1) a constant initial condition for a single drug delivery experiment; and (2) a recursive delivery process where the liquid medium is replaced with fresh liquid after a fixed period of time, causing a stepped delivery rate.

2.1. Single Drug Delivery Experiments

Consider a unidirectional diffusion of the solute in the solid matrix, described by the following partial differential equation:
C x , t t = D 2 C x , t x 2 , t > 0 0 < x < L
where D is the diffusivity and L is the length. The initial condition is
C x , 0 = C i n i
and boundary conditions are
C 0 , t x = 0
D C L , t x = k C C L , t γ C l t
where C l t is the solute concentration in the liquid, k C is the mass transfer coefficient in the liquid, and γ is the solute solid–liquid partition coefficient associated with the equilibrium. Condition (2) indicates that the initial concentration in the solid matrix is homogeneous and equal to C i n i . Condition (3) appears due to the symmetry of diffusion in the middle of the solid, while condition (4) represents the solid–liquid interface. The accumulation of the solute in the liquid can be easily obtained from the mass balance,
V l d C l t d t = A k C C L , t γ C l t
where A is the area of mass transport and V l is the volume of the liquid. The total solute mass in the solid matrix and the liquid is constant and equal to
m i n i = A 0 L C x , t d x + V l C l t
where m i n i is the initial solute mass in the solid matrix and, assuming that initially the solute concentration in the liquid is zero C l 0 = 0 , it is equivalent to
m i n i = A 0 L C x , 0 d x = A L C i n i .
Notice that the time derivative of Equation (6) provides
V l d C l t d t = A 0 L C x , t t d x
and considering Equation (1) and boundary condition (3), the previous equation becomes
V l d C l t d t = D A C L , t x ,
which agrees with the substitution of Equation (4) in Equation (5). The mass of the solute in the liquid is therefore
C l t = m i n i m t V l = A L C i n i 0 L C x , t d x V l
where m t = A 0 L C x , t d x is the total mass of the solute in the solid matrix. Finally, replacing Equation (8) in boundary condition (4), it follows that
L Sh C L , t x + C L , t = β m i n i A L 1 L 0 L C x , t d x
where β = γ A L V l and Sh = k C D / L are the Sherwood number, which represents the ratio of convective mass transfer to the rate of diffusive mass transport.

2.2. Recursive Delivery Experiment

Now, let us consider an iterative process where for each fixed period of time δ the liquid is replaced by fresh liquid without solute. Then, the iterative process can be modeled by the diffusion equation,
C k x , t t = D 2 C k x , t x 2 , k δ < t k + 1 δ 0 < x < L
with an initial condition that depends on the concentration of the previous iteration:
C k x , k δ = C k 1 x , k δ
and the following boundary conditions:
C k 0 , t x = 0
L Sh C L , t x + C k L , t = β m 0 , k A L 0 L C k x , t d x L k δ < t k + 1 δ
where k = 1 , 2 , 3 , and the initial mass of the solute in the solid at iteration k depends on the final concentration of iteration k 1 , i.e.,
m 0 , k = A 0 L C k 1 x , k δ d x .
To analyze the iterative behavior, it is proposed that a shift in time τ k = t k δ be defined in such a way that the previous model becomes
C k x , τ k τ k = D 2 C k x , τ k x 2 , τ k > 0 0 < x < L
C k x , 0 = C k 1 x , τ k 1 τ k 1 = δ
C k 0 , τ k x = 0
L Sh C k L , τ k x + C k L , τ k = β m 0 , k A L 0 L C k x , τ k d x L
for k = 1 , 2 , 3 , . Then, the analytical solution of the models for the two scenarios is presented in the following section using Laplace transformation [6].

3. Analytical Solutions

3.1. Single Drug Delivery Experiment

To solve the problem described by Equations (1)–(3) and (9), let us define the Laplace transformation:
L C x , t = 0 C x , t e s t d t = C ^ x , s
therefore, the Laplace transformation of Equation (1) is
s C ^ x , s C i n i = D d 2 C ^ x , s d x 2 , 0 < x < L
while the Laplace transformations of the boundary conditions (Equations (3) and (9)) are
d C ^ 0 , s d x = 0 ,
L Sh d C ^ L , s d x + C ^ L , s = β m i n i A L s 1 L 0 L C ^ x , s d x .
The solution of Equation (19) together with Conditions (20) and (21) is
C ^ x , s = C i n i s 1 ρ cosh x s D cosh L s D + 1 L D s L 2 s D Sh + β sinh L s D
where ρ = 1 + β 1 m i n i C i n i A L . Notice that, for this particular case ρ = 1 , because m i n i = A L C i n i ; however, here, we use ρ , as the general solution for ρ 1 is used for the recursive problem analyzed latter.
The concentration profile can be obtained by applying the inverse Laplace transformation to Equation (22). The inverse of the first term can be obtained with the formula L 1 1 s = 1 , while for the second term a more elaborated procedure is required.
We define
p s = cosh x L L 2 s D
and
q s = cosh L 2 s D + D L 2 s β + 1 Sh L 2 s D sinh L 2 s D .
where for s n = D α n L 2 , it holds that
q s n = cos α n + β α n α n Sh sin α n .
Thus, the values of α n at which q s n = 0 are provided by the characteristic equation
tan α n = α n α n 2 Sh β .
according to the following formula:
L 1 p s q s = n = 1 p s n q s n exp s n t
where p s and q s are infinite series of s, q s is the derivative of q s with respect to s, and s n are the infinite values of s satisfying the equation q s = 0 ; thus, the following inverse transformation is obtained
L 1 p s q s = 2 D L 2 n = 1 α n 3 cos α n x L exp D α n L 2 t β + α n 2 1 + 1 Sh + β α n 2 Sh 2 sin α n .
Finally, for any arbitrary function f s = L f t , the following inverse Laplace transformation holds: L 1 f s s = 0 t F τ d τ ; then,
L 1 1 s · p s q s = 2 n = 1 cos α n x L 1 exp D α n L 2 t α n 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β α n sin α n
and the concentration becomes the following Fourier series [7]:
C x , t = C i n i 1 2 ρ n = 1 cos α n x L 1 exp D α n L 2 t α n 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β α n sin α n .
The relationship between the mass at any time and the initial mass of solute in the solid matrix is
m t m i n i = 1 L 0 L C x , t C i n i d x
and considering the solution, this integral is
0 L C x , t C i n i d x L = 1 2 ρ n = 1 1 exp D α n L 2 t α n 2 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β
while the amount of mass of solute that has left the solid matrix and diffused to the liquid is
1 m t m i n i = 2 ρ n = 1 1 exp D α n L 2 t α n 2 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β
At steady state, the concentration in the solid matrix does not depend on time; thus, C x , t becomes C * x . Therefore, Equation (1) reduces to
D d 2 C * x d x 2 = 0
while conditions (3) and (4) become
d C * 0 d x = 0
L Sh d C * L d x + C * L = β m i n i A L 1 L 0 L C * x d x
The solution of Equation (26) together with the boundary conditions (27) and (28) is
C * x = β 1 + β m i n i A L = γ C l ,
where C l , is the steady state concentration of the solute in the liquid. Notice that the steady state concentration profile of the solute in the solid matrix is constant.
Now, considering the transient solution (25) and the solution at steady state (29), given that C * x = lim t C x , t it is straightforward to verify that the following Fourier series can be expressed as
2 n = 1 cos α n x L α n 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β α n sin α n = 1 1 + β
and solution (25) is equivalent to
C x , t = β 1 + β m i n i A L + 2 ρ C i n i n = 1 cos α n x L exp D α n L 2 t α n 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β α n sin α n .
The integral of this concentration is
m t m i n i = 1 L 0 L C x , t C 0 d x = β 1 + β m i n i A L C i n i + 2 ρ n = 1 exp D α n L 2 t α n 2 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β ,
and the mass of the solute in the liquid is
m l t m i n i : = 1 m t m i n i = ρ 1 1 + β 2 n = 1 exp D α n L 2 t α n 2 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β .

3.2. Recursive Delivery Experiment

In order to solve this problem, we first state that
C k x , τ k = C k 1 x , τ k 1 + δ + W k x , τ k ,
where W k x , τ k is an auxiliary variable and C k 1 x , τ k 1 + δ is the solution at the previous iteration. Therefore, the substitution of Equation (33) in Equation (14) and conditions (15)–(17) produces the following problem for W k x , τ k :
W k x , τ k τ k = D 2 W k x , τ k x 2 , τ k > 0 0 < x < L
W k x , 0 = 0
W k 0 , τ k x = 0
L Sh W k L , τ k x + W k L , τ k + β 0 L W k x , τ k d x L = β m 0 , k m 0 , k 1 A L
where
m 0 , k = A 0 L C k 1 x , τ k 1 τ k 1 = δ d x , m 0 , k 1 = A 0 L C k 1 x , τ k 1 τ k 1 = 0 d x .
Notice that Equations (34)–(37) are analogue to Equations (1)–(3) and (9) with C i n i = 0 and m i n i = m 0 , k m 0 , k 1 ; therefore, the solution is
W k x , τ k = 2 β m 0 , k m 0 , k 1 A L n = 1 cos α n x L 1 exp D α n L 2 τ k α n 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β α n sin α n
while considering the identity in (30), this function becomes
W k x , τ k = β m 0 , k m 0 , k 1 A L 1 W 0 x , τ k
where
W 0 x , t = β 1 + β + 2 n = 1 cos α n x L exp D α n L 2 t α n 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β α n sin α n .
The concentration provided in Equation (33) becomes
C k x , τ k = C k 1 x , τ k 1 + δ + β m 0 , k m 0 , k 1 A L 1 W 0 x , τ k ,
which in terms of the original time variable is
C k x , t = C k 1 x , t + β m 0 , k m 0 , k 1 A L 1 W 0 x , t k δ , k δ < t k + 1 δ 0 < x < L .
Notice that for k = 0 , the solution is obtained through Equations (25) or (31)
C 0 x , t = C i n i W 0 x , t , 0 < t δ 0 < x < L ;
for k = 1 , it holds that
C 1 x , t = C i n i W 0 x , t + β m 0 , 1 A L C i n i 1 W 0 x , t δ , δ < t 2 δ 0 < x < L ,
while for k = 2 ,
C 2 x , t = C i n i W 0 x , t + β m 0 , 1 A L C i n i 1 W 0 x , t δ + β m 0 , 2 m 0 , 1 A L 1 W 0 x , t 2 δ , 2 δ < t 3 δ 0 < x < L .
and in general,
C k x , t C i n i = W 0 x , t + β i = 1 k m 0 , i m 0 , 0 m 0 , i 1 m 0 , 0 1 W 0 x , t i δ , k δ < t k + 1 δ 0 < x < L
where m 0 , 0 = A L C i n i = m i n i .
The total fraction of solute mass inside the solid matrix with respect to the initial solute mass at each instant is
m k t m 0 , 0 = A 0 L C k x , t d x m 0 , 0 = M 0 t + β i = 1 k m 0 , i m 0 , 0 m 0 , i 1 m 0 , 0 1 M 0 t i δ
where
M 0 t = 1 L 0 L W 0 x , t d x = β 1 + β + 2 n = 1 exp D α n L 2 t α n 2 1 + 1 2 β + α n 2 / Sh Sh + 1 + β β
therefore,
m 0 , 1 m 0 , 0 = M 0 δ
and
m 0 , k m 0 , 0 = m k 1 k δ m 0 , 0 = M 0 k δ + β i = 1 k 1 m 0 , i m 0 , 0 m 0 , i 1 m 0 , 0 1 M 0 k i δ , k = 2 , 3 ,
Finally, the mass of the solute that has been delivered from the solid matrix is
m l , k t m 0 , 0 : = 1 m k t m 0 , 0 = 1 M 0 t β i = 1 k m 0 , i m 0 , 0 m 0 , i 1 m 0 , 0 1 M 0 t i δ .

4. Case Study

4.1. Materials and Methods

Surface sections ( 1 × 1 cm 2 ) from a roughness breast implant purchased from Mentor® brand (USA) were polymerized with (2-Hydroxypropyl)-b-cyclodextrin and citric acid according to the protocol of [8]. Samples were loaded with a concentrated solution of Rose Bengal (RB) stain for 12 h. The RB release profile, C l t , was obtained as a mean of triplicate experiments. Loaded samples were placed into vials filled with phosphate-buffered saline (PBS) at 37 °C in a horizontal shaker (100 rpm) (NB-2005LN Biotek, Winooski, VT, USA). In order to maintain perfect sink conditions, the PBS medium of the samples was changed every 24 h, for ten times in total. The RB content in the withdrawn bulk PBS was analyzed by UV-spectrophotometry (Evolution 220 model, Thermo Scientific, Waltham, MA, USA) at 545 nm [9]. Reactive chemicals were obtained from Sigma-Aldrich® (St. Louis, MO, USA).

4.2. Model Prediction

Figure 2a shows the RB release profile for the first 24 h of the experiment as m l t / m i n i vs. t. As a first approach, these experimental data were used to estimate the three parameters of the model presented in Section 2.1  D , β , Sh , with a least-squares method using a derivative-free algorithm [10] to minimize the sum of the squares of the error between the data and the prediction of Equation (32). As a result, a very large value was obtained for the Sherwood number, Sh = k C D / L 10 5 , suggesting that diffusion controls the interface transport for the experimental conditions; therefore, the characteristic Equation (23), the concentration profile (31), and the release profile (32) can be simplified to tan α n = α n / β :
C x , t C i n i = β 1 + β + 2 n = 1 cos α n x L exp D α n L 2 t α n + 1 + β β α n sin α n
and
m l t m i n i = 1 1 + β 2 n = 1 exp α n 2 D L 2 t α n 2 + 1 + β β ,
respectively. This simplified model contains only two parameters, D / L 2 , β , the numerical values of which, D / L 2 = 4.102 × 10 2 h 1 and β = 6.095 × 10 2 , were estimated using the same least-square method [10]. Finally, the same parameters were used to predict the behavior of the recursive delivery experiment, where the PBS medium of the samples was changed every 24 h for ten days (see Figure 3a).

5. Discussion

The proposed models for single and recursive release experiments presented in Section 2.1 and Section 2.2, respectively, are general and allow us to clarify the predominant step in the control of drug release, This is because these models consider the drug diffusion in the solid matrix, the convective phenomenon from the solid matrix to the liquid, where the drug is delivered, and the solid–liquid drug distribution equilibrium. In particular, considering the experimental conditions of the study case presented in Section 4, diffusion dominates over convection thanks to the mixing conditions established by the horizontal shaker (with 100 rpm) in our experiments. For this reason, we was obtained Sh 1 0 , reducing the number of parameters to two, namely, D / L 2 , β ; however, convection might be relevant as well for in situ applications.
In the Higuchi model [3], the characteristic values are of the form n 1 2 π , leading to the approximation M t / M = k t 1 / 2 . However, the assumption of solid–liquid interface equilibrium leads to the characteristic Equation (23), which depends on both Sh and β and the characteristic values of which must be numerically computed. In particular, for Sh 1 = 0 and β = 6.095 × 10 2 the first characteristic values are α 1 0.6 π , α 2 1.54 π , α 3 2.524 π , α 4 3.518 π ,…, which tends towards α n n 1 2 π as n increases (see Figure 4). The deviation from n 1 2 π of the first characteristic values produces a response that is non necessarily proportional to t 1 / 2 , as in the Higuchi model [3].
The experimental evidence shown in Figure 2a and Figure 3a supports the use of the solid–liquid interface equilibrium assumption, C e q = γ C l , because at the end of the first 24 h the concentration in the liquid seems to reach a constant value (see Figure 2a); it was only after the PBS medium was exchanged for a fresh one that the drug release continued (see Figure 2a). Condition (9) or its equivalent for the recursive delivery (17) includes the solid–liquid interface equilibrium assumption, and these conditions indicate that the rate of drug delivery to the liquid is time-varyiant, as it depends on the total mass of the drug in the liquid at the same instant. This concentration is proportional to the total mass inside the solid matrix due to the total mass balance of the total solid–liquid system. As a result of these conditions, the delivered mass with respect to the initial mass within the solid matrix tends to 1 + β 1 as time tends to infinity (see Equation (32)), which is consistent with the experimental data, i.e., m l / m i n i 0.621 for β = 6.095 × 10 2 . Thus, for the particular experimental study case, when the drug concentration in the liquid is low, diffusion controls the delivery rate; however, as this concentration increases, the equilibrium plays a greater role, becoming the controlling mechanism as time increases.
The proposed models for single or recursive delivery are able to describe the delivered mass profile (see Figure 2a and Figure 3a), that is, the usual experimental measurement, and allows prediction of the concentration profile within the solid matrix, as shown in Figure 2b and Figure 3b. Notice that for the particular experimental study case, after each 24 h the profile concentration in the solid matrix is approximately constant, because the solid–liquid interface equilibrium has been reached. In addition, in order to validate the analytical solutions provided in Equation (32) for single release experiments and Equation (40) for recursive release experiments, numerical solutions of the proposed models were determined using an orthogonal collocation method; Figure 2a and Figure 3a show that the numerical and analytical solutions overlap.

6. Conclusions

In this work, a mathematical model to describe drug delivery from polymer coatings on implants was proposed. The model considers a unidirectional recursive diffusion process which follows Fick’s second law and considers the convective phenomena from the polymer matrix to the liquid where the drug is delivered as well as the polymer–liquid drug distribution equilibrium. The resulting model is solved using Laplace transformation for two scenarios: (1) a constant initial condition for a single drug delivery experiment; and (2) a recursive delivery process where the liquid medium is replaced with fresh liquid after a fixed period of time, causing a stepped delivery rate. Finally, the study case shows that these models can satisfactorily reproduce the experimental data, clarifying the predominant step in the control of drug release.
The proposed models for single and recursive release experiments presented in Section 2.1 and Section 2.2 have linear partial differential equations and boundary conditions; therefore, it was possible to obtain the analytical solutions provided in Equations (25), (31) and (32) for single release experiments, and Equations (39) and (40) for recursive release experiments. However, numerical methods are another option to solve these problems, and they become mandatory when extra phenomena need to be included in the drug release models in order to obtain satisfactory predictions, for instance, when swelling and degradation of the solid matrix become relevant, causing moving boundaries, or when a drug has more complex interactions such as biodegradation or adsorption, which produce non-linear terms in the model due to the reaction kinetics or adsorption isotherms.

Supplementary Materials

The following supporting information can be downloaded at: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/math10132171/s1.

Author Contributions

Conceptualization, J.P.G.-S. and J.S.-A.; methodology, J.H.-M.; software, J.P.G.-S. and J.S.-A.; validation, J.P.G.-S., J.H.-M. and R.H.-M.; formal analysis, J.P.G.-S. and J.S.-A.; investigation, J.S.-A.; resources, J.H.-M.; data curation, J.S.-A.; writing—original draft preparation, J.P.G.-S., J.H.-M. and R.H.-M.; writing—review and editing, J.H.-M.; visualization, R.H.-M.; supervision, J.H.-M.; project administration, J.H.-M.; funding acquisition, J.H.-M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by FONDECYT, Chile, grant number 11180395.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in Supplementary Materials.

Acknowledgments

J.S.-A. acknowledges support by ANID Master’s degree scholarship 22212018 and Effective Access to Higher Education Programme, PACE-UFRO. J.P.G.-S. and R.H.-M. acknowledge the grants (SNI-43658 and SNI-149478) provided by the Mexican National Council of Science and Technology (CONACyT) through the program Sistema Nacional de Investigadores (SNI).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Quarterman, J.C.; Geary, S.M.; Salem, A.K. Evolution of drug-eluting biomedical implants for sustained drug delivery. Eur. J. Pharm. Biopharm. 2021, 159, 21–35. [Google Scholar] [CrossRef] [PubMed]
  2. Siepmann, J.; Siepmann, F. Mathematical modeling of drug delivery. Int. J. Pharm. 2008, 364, 328–343. [Google Scholar] [CrossRef] [PubMed]
  3. Higuchi, T. Mechanism of sustained-action medication. Theoretical analysis of rate of release of solid drugs dispersed in solid matrices. J. Pharm. Sci. 1963, 52, 1145–1149. [Google Scholar] [CrossRef]
  4. Ritger, P.L.; Peppas, N.A. A simple equation for description of solute release I. Fickian and non-fickian release from non-swellable devices in the form of slabs, spheres, cylinders or discs. J. Control. Release 1987, 5, 23–36. [Google Scholar] [CrossRef]
  5. Fu, Y.; Kao, W.J. Drug release kinetics and transport mechanisms of non-degradable and degradable polymeric delivery systems. Expert Opin. Drug Deliv. 2010, 7, 429–444. [Google Scholar] [CrossRef]
  6. Churchill, R.V. Operational Mathematics, 3rd ed.; McGraw-Hill: New York, NY, USA, 1972. [Google Scholar]
  7. Brown, J.W.; Churchill, R.V. Fourier Series and Boundary Value Problems, 8th ed.; McGrawHill: New York, NY, USA, 2012. [Google Scholar]
  8. Hernandez-Montelongo, J.; Naveas, N.; Degoutin, S.; Tabary, N.; Chai, F.; Spampinato, V.; Ceccone, G.; Rossi, F.; Torres-Costa, V.; Manso-Silvan, M.; et al. Porous silicon-cyclodextrin based polymer composites for drug delivery applications. Carbohydr. Polym. 2014, 110, 238–252. [Google Scholar] [CrossRef]
  9. França, C.G.; Plaza, T.; Naveas, N.; Santana, M.H.A.; Manso-Silván, M.; Recio, G.; Hernandez-Montelongo, J. Nanoporous silicon microparticles embedded into oxidized hyaluronic acid/adipic acid dihydrazide hydrogel for enhanced controlled drug delivery. Microporous Mesoporous Mater. 2021, 310, 110634. [Google Scholar] [CrossRef]
  10. Lagarias, J.C.; Reeds, J.A.; Wright, M.H.; Wright, P.E. Convergence Properties of the Nelder–Mead Simplex Method in Low Dimensions. SIAM J. Optim. 1998, 9, 112–147. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Scheme of the drug delivery system.
Figure 1. Scheme of the drug delivery system.
Mathematics 10 02171 g001
Figure 2. Single drug delivery. (a) Comparison of the experimental data and prediction for the delivered mass using Equation (32). (b) Prediction of the drug concentration profile in the solid matrix using Equation (31). The parameters are D / L 2 = 4.102 × 10 2 h 1 , β = 6.095 × 10 2 , and Sh 1 = 0 .
Figure 2. Single drug delivery. (a) Comparison of the experimental data and prediction for the delivered mass using Equation (32). (b) Prediction of the drug concentration profile in the solid matrix using Equation (31). The parameters are D / L 2 = 4.102 × 10 2 h 1 , β = 6.095 × 10 2 , and Sh 1 = 0 .
Mathematics 10 02171 g002
Figure 3. Recursive drug delivery. (a) Comparison of the experimental data and prediction for the delivered mass using Equation (40). (b) Prediction of the drug concentration profile in the solid matrix using Equation (39). The parameters are D / L 2 = 4.102 × 10 2 h 1 , β = 6.095 × 10 2 , and Sh 1 = 0 .
Figure 3. Recursive drug delivery. (a) Comparison of the experimental data and prediction for the delivered mass using Equation (40). (b) Prediction of the drug concentration profile in the solid matrix using Equation (39). The parameters are D / L 2 = 4.102 × 10 2 h 1 , β = 6.095 × 10 2 , and Sh 1 = 0 .
Mathematics 10 02171 g003
Figure 4. Characteristic values, α n , for β = 6.095 × 10 2 .
Figure 4. Characteristic values, α n , for β = 6.095 × 10 2 .
Mathematics 10 02171 g004
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hernandez-Montelongo, R.; Salazar-Araya, J.; Hernandez-Montelongo, J.; Garcia-Sandoval, J.P. Mathematical Modeling of Recursive Drug Delivery with Diffusion, Equilibrium, and Convection Coupling. Mathematics 2022, 10, 2171. https://0-doi-org.brum.beds.ac.uk/10.3390/math10132171

AMA Style

Hernandez-Montelongo R, Salazar-Araya J, Hernandez-Montelongo J, Garcia-Sandoval JP. Mathematical Modeling of Recursive Drug Delivery with Diffusion, Equilibrium, and Convection Coupling. Mathematics. 2022; 10(13):2171. https://0-doi-org.brum.beds.ac.uk/10.3390/math10132171

Chicago/Turabian Style

Hernandez-Montelongo, Rosaura, Javiera Salazar-Araya, Jacobo Hernandez-Montelongo, and Juan Paulo Garcia-Sandoval. 2022. "Mathematical Modeling of Recursive Drug Delivery with Diffusion, Equilibrium, and Convection Coupling" Mathematics 10, no. 13: 2171. https://0-doi-org.brum.beds.ac.uk/10.3390/math10132171

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop