Abstract
A theoretical model of the drug release process from polymeric microparticles (a particular type of polymer matrix), through dispersive fractal approximation of motion, is built. As a result, the drug release process takes place through cnoidal oscillations modes of a normalized concentration field. This indicates that, in the case of long-time-scale evolutions, the drug particles assemble in a lattice of nonlinear oscillators occur macroscopically, through variations of drug concentration. The model is validated by experimental results.
1. Introduction
Polymer matrices can be produced in one of the following forms: micro/nanoparticles, micro/nanocapsules, hydrogels, films, and patches. Due to the multitude of biocompatible polymers in the experimental protocol, drug proper delivery via many administration routes occurs. No matter what their form might be, drug carrier polymer matrices should have the following characteristics: biocompatibility, biodegradability, and controlled release capacity. The last one refers to the relationship between the efficient, nontoxic drug administration and therapeutic window type concentration, that is, minimum concentration is required to produce the wanted effect, but in the case of high levels, a toxic barrier occurs.
Given the importance of the released drug concentration, numerous studies have been performed with the purpose of identifying the mathematical function that describes time dependence. Many papers show how various factors, such as polymer molecular weight [1, 2], polymer chemistry, monomer ratios [3, 4], pH of release media, additives to the release media [5, 6], and particle size [7], affect the release kinetics. At the same time, certain phenomena appearing in the release process have been studied. Of these, we mention (in approximate order of their occurrence) polymer swelling and degradation [8–11], drug dissolution and diffusion [12, 13], and above all, permanent chemical and physical interaction among components (drug, polymer, and release medium). Since all these phenomena are not independent, their analysis becomes complicated; consequently, it will not be possible to treat them separately and cumulate the effects. For example, microparticle morphology changes due to polymer degradation, their surfaces becoming highly porous. This will lead to increased diffusion coefficients and hence certain connected phenomena, such as polymer degradation and drug diffusion [7].
The multitude of phenomena and dependencies occurring in drug release process as well as numerous structural entities (polymer, drug, and release medium) will turn the system into a complex one. Consequently, the complete theoretical analysis becomes more difficult in terms of performing.
Nevertheless, significant amount of work has been accomplished in mathematical modeling, with the purpose of predicting the concentration of the released drug and providing the analysis of fundamental processes that govern release. Higuchi [14] was among the first who produced a drug release model from nonswelling and nondissolving polymer matrices, assuming that such phenomenon is purely controlled through diffusion. A number of other models have also been proposed in order to predict drug release in the case of erosion [9, 11], swelling [8], and dissolution [13] influenced processes. These mathematical patterns have chosen only two phenomena, with the purpose of simplifying mathematical modeling, which, otherwise, proves to be quite difficult.
That is the main reason why it is necessary to use alternative approaches with reduced number of the approximations. One of such possible approaches is the fractal one [15, 16]. Its use is justified by natural and synthetic polymers that have been included in the category of fractional-dimensioned objects whose structures and behaviour can be described by means of fractal geometry [17, 18]. Moreover, it has been observed that the dynamics of drug release systems is a fractal one, because, in spite of complex phenomena and factors, mathematical expressions describing drug release kinetics from a variety of polymer matrices are power type laws (Higuchi [14] for nonswelling and nondissolving polymer, Ritger and Peppas [19] for nonswellable polymer in the form of slabs, spheres, cylinders, or discs, Peppas Sahlin [20] for solute release, Alfrey et al. [21] for diffusion in glassy polymers, etc.) specific for the fractal system evolution [22]. At the same time, it is quite important to emphasize that correlation of experimental data with the above-mentioned laws revealed good correlation in the first part (approximate 60%) of the release kinetics, the correlation coefficient decreasing according to time evolution.
Studies on the release from different types of systems (HPMC matrices [23], inert porous matrices [24], and sponges [25, 26]) have been performed. Such approach analyzes drug release kinetics through Monte Carlo simulation. In this perspective, release systems are considered as three-dimensional lattices with leak sites located at the boundaries of the lattice pattern. Particles are free to move inside the porous network according to the random walk model of the Fickian diffusion (the moving particles act as hard spheres colliding with each other and having no possibility to mutually penetrate).
The first studies started with simplifying approximations. Kosmidis et al. [23] consider that porosity has a constant value. Later on, Villalobos et al. [24] improved the model, assumed that network porosity behaves dynamically, and considered the effects of drug spatial distribution and initial drug concentration. All these approaches proved the validity of Weibull function (a continuous probability distribution function) for the entire release kinetics and consequently eliminate Peppas’ temporal limitation of the equation and criticism lacking kinetic basis and physical nature of parameters [27].
Our new approach considers the entire system (drug-loaded polymer matrix in the release environment) as a type of “fluid” totally lacking interaction or neglecting physical interactions among particles. At the same time, the induced complexity is replaced by fractality. This will lead to particles moving on certain trajectories called geodesics within fractal space. This assumption represents the basis of the fractal approximation of motion in scale relativity theory (SRT) [28, 29], leading to a generalized fractal “diffusion” equation that can be analyzed in terms of two approximations (dissipative and dispersive).
The comparison between dissipative approximation (with dominant convective and dissipative processes) and the dispersive one allows theoretical demonstration of Weibull function that best describes the behaviour of drug release systems at short time scales. These phenomena will be the object of subsequent analysis since they are responsible for certain types of behaviour and characterized by high degree of nonlinearity in drug release systems.
This paper is structured as follows: theoretical model (Section 2), experimental results that validate the theoretical model (Section 3), and conclusions (Section 4).
2. Theoretical Model
2.1. Consequences of Nondifferentiability
We suppose that the drug release process takes place on continuous, but nondifferentiable curves (fractal curves). Then, nondifferentiability implies [28–30] the following.(i)A continuous and a nondifferentiable curve (or almost nowhere differentiable) is explicitly scale dependent, and its length tends to infinity, when the scale interval tends to zero. In other words, a continuous and nondifferentiable space is fractal, and in the general meaning Mandelbrot used this concept [15];(ii)Physical quantities will be expressed through fractal functions, namely, through functions that are dependent both on coordinate field and resolution scale. The invariance of the physical quantities in relation with the resolution scale generates special types of transformations, called resolution-scale transformations. In what follows, we will explain the above statement.
Let be a fractal function in the interval , and let the sequence of values for be
Let us note by the broken line that connects the points
We can now say that is a -scale approximation.
Let us now consider as a -scale approximation of the same function. Since is everywhere almost self-similar, if and are sufficiently small, both approximations and must lead to the same results; in the particular case, a fractal phenomenon is studied through approximation. By comparing the two cases, one notices that scale expansion is related to the increase of , according to an increase of (see Figure 1). But, in this case, we have a situation in which we can consider the infinitesimal-scale transformation as being
Such transformation in the case of function leads to respectively, if we limit ourselves to a first-order approximation:
Moreover, let us notice that for an arbitrary but fixed , we obtain a situation in which (2.6) can be written as follows:
Therefore, we can introduce the dilatation operator
At the same time, relation (2.9) shows that the intrinsic variable of resolution is not , but .
The fractal function is explicitly dependant on the resolution ; therefore, we have to solve the differential equation where is now an unknown function. The simplest explicit suggested form for is linear dependence [29] in which case the differential equation (2.10) takes the form
Hence, by integration and substitutingwe obtain
This solution is independent as compared to parameterization on fractal curve.
We can now generalize the previous result by considering that is dependent on parameterization of the fractal curve. If characterizes the position on the fractal curve, then, following the same algorithm as above, the solution will be as a sum of two terms, that is, both classical and differentiable (depending only on position) and fractal, nondifferentiable (depending on position and, divergently, on ) where is a function depending on parameterization of the fractal curve.
The following particular cases are to be considered.
In asymptotic small-scale regime, is constant (with no scale dependence) and power-law dependence on resolution is obtained:
At this stage, some power laws should also be considered, namely, those equations describing drug release kinetics from a different type of polymer matrix [14, 19–21]. Consequently, through the appropriate correspondence among quantities from (2.16a) and (2.16b) and those from drug release processes, we will obtain the following: (a) Higuchi law: where and are cumulative amounts of drug release at time and infinity, respectively, and is a constant characteristic of the system [14];(b)Peppas law: where is an experimentally obtained parameter, and is a real number geometrically related to the system and to drug release mechanism. The value is used to characterize different release mechanisms, that is, indicates a Fickian diffusion. In their turn, different from 0.5 values refer to mass transport according to non-Fickian model [31]. This equation is a generalization of a square-root time law and an approximation for short times of Weibull function.
In these expressions, we recognize the standard form of self-similar fractal behaviour with fractal dimension , which has already been used for accurately describing many physical and biological systems [15]. The topological dimensions are here , since we deal with length, but this can be easily generalized to surfaces () and volumes (). Therefore, such result is not a consequence of postulation or deduction, but an aftermath of first principle theoretical analysis.
Considering that the resolution is a length, , the scale-dependent length is given, by definition, by the law where is a scale characteristic length, and the exponent is identified with .
Now, in the above solution, one may use time as parameter, and if one constantly moves along the curve, one obtains . Then, a differential version of the above relation will be so that the following fundamental relation among space and time elements on a fractal curve or function is obtained: In other words, they are differential elements of different orders.
In the asymptotic big-scale regime , is constant (with no scale dependence), and, in terms of resolution, one obtains an independent law
Particularly, if are the coordinates in given space, we can write
In this situation, becomes a highly fluctuating function which can be described by stochastic process, while represents (according to previous description) the difference between fractal and topological dimensions. The result is a sum of two terms, a classical, differentiable one (dependent only on the position) and a fractal, nondifferentiable one (dependent both on the position and, divergently, on ). This represents the importance of the above analysis.
By differentiating these two parts, we obtain where is the classical differential element, and is a differential fractal one.
(iii) There is infinity of fractal curves (geodesics) relating to any couple of points (or starting from any point) and applied for any scale. The phenomenon can be easily understood at the level of fractal surfaces, which, in their turn, can be described in terms of fractal distribution of conic points of positive and negative infinite curvature. As a consequence, we have replaced velocity on a particular geodesic by fractal velocity field of the whole infinite ensemble of geodesics. This representation is similar to that of fluid mechanics [32] where the motion of the fluid is described in terms of its velocity field , density , and, possibly, its pressure. We will, indeed, recover the fundamental equations of fluid mechanics (Euler and continuity equations), but we will write them in terms of a density of probability (as defined by the set of geodesics) instead of a density of matter and adding an additional term of quantum pressure (the expression of fractal geometry).
(iv) The local differential time invariance is broken, so the time derivative of the fractal field can be written as twofold:
Both definitions are equivalent in the differentiable case . In the nondifferentiable situation, these definitions are no longer valid, since limits are not defined anymore. Fractal theory defines physics in relationship with the function behavior during the “zoom” operation on the time resolution , here identified with the differential element (substitution principle), which is considered an independent variable. The standard field is therefore replaced by fractal field , explicitly dependent on time resolution interval, whose derivative is not defined at the unnoticeable limit . As a consequence, this leads to the two derivatives of the fractal field as explicit functions of the two variables and ,
Notation “+” corresponds to the forward process, while “−” to the backward one.
(v) Let be a point of the fractal curve, and let us consider a line which starts from this point. Let be the first intersection of this line with the fractal curve. By , we denote the components of the vector PM, to the right of the line (d), and by the components of the vector PM′, to the left of the line (d)—see Figure 2.
If we consider all the lines (segments) starting from , we denote the average of these vectors by , that is,
Since, according to (2.24), we can write
it results that
(vi) The differential fractal part satisfies, according to (2.21), the fractal equation
where and are some constant coefficients, and is a constant fractal dimension. We note that the use of any Kolmogorov or Hausdorff [15, 28, 33–35] definitions can be accepted for fractal dimension, but once a certain definition is admitted, it should be used until the end of analyzed dynamics.
(vii) The local differential time reflection invariance is recovered by combining the two derivatives, and , in the complex operator
Applying this operator to the “position vector,” a complex speed yields with
The real part, , of the complex speed represents the standard classical speed, which does not depend on resolution, while the imaginary part, , is a new quantity coming from resolution-dependant fractal.
2.2. Covariant Total Derivative in Drug Release Mechanism
Let us now assume that curves describing drug release (continuous but nondifferentiable) are immersed in a 3-dimensional space, and that of components is the position vector of a point on the curve. Let us also consider the fractal field and expand its total differential up to the third orderwhere only the first three terms were used in Nottale’s theory (i.e., second-order terms in the motion equation). Relations (2.34a) and (2.34b) are valid in any point both for the spatial manifold and for the points on the fractal curve (selected in relations (2.34a) and (2.34b)). Hence, the forward and backward average values of these relations take the form
The following aspects should be mentioned: the mean value of function and its derivatives coincide with themselves, and the differentials and are independent; therefore, the average of their products coincides with the product of averages. Consequently, (2.35) and (2.36) become
or more, using (2.28a) and (2.28b) with characteristics (2.29a) and (2.29b),
Even if the average value of the fractal coordinate is null (see (2.29a) and (2.29b)), for higher order of fractal coordinate average, the situation can still be different. Firstly, let us focus on the averages and . If , these averages are zero due to the independence of and . So, using (2.30a) and (2.30b), we can write
Then, let us consider the averages and . If , these averages are zero due to independence of on and . Now, using (2.30a) and (2.30b), we can write
Then, (2.38a) and (2.38b) may be written as follows:
If we divide by and neglect the terms containing differential factors (for details on the method, see [36, 37]), (2.41a) and (2.41b) are reduced to
These relations also allow us to define the operator
Under these circumstances, let us calculate . Taking into account (2.43a), (2.43b), (2.31), and (2.32), we will obtain
This relation also allows us to define the fractal operator
Particularly, by choosingthe fractal operator (2.45) takes the usual form
We now apply the principle of scale covariance and postulate that the passage from classical (differentiable) to “fractal” mechanics can be implemented by replacing the standard time derivative operator, , with the complex operator (this results in a generalization of Nottale’s [28, 29] principle of scale covariance). Consequently, we are now able to write the diffusion equation in its covariant form
This means that at any point on a fractal path, the local temporal , the nonlinear (convective), , the dissipative, , and the dispersive, , terms keep their balance.
The dissipative approximation was applied for the drug release processes, and the result was a Weibull type function that was analyzed in [38, 39]. In what follows, we will focus on dispersive approximation.
2.3. The Dispersive Approximation
Let us now consider that, in comparison with dissipative processes, convective and dispersive processes are dominant ones. Consequently, we are now able to write the diffusion equation in its covariant form, as a Korteweg de Vries type equation
If we separate the real and imaginary parts from (2.49), we will obtain
By adding them, the fractal diffusion equation is
From (2.50b), we see that, at fractal scale, there will be no field gradient.
Assuming that with (in systems with self-structuring processes, the speed fluctuations induced by fractal/nonfractal are proportional with the concentration field [22]), in the particular one-dimensional case, (2.51) will lack parameters
and normalizing conditions take the form
In relations (2.52a), (2.52b), (2.52c), and (2.53), corresponds to a characteristic pulsation, to the inverse of a characteristic length, and to balanced concentration.
Through substitutions
(2.54), by double integration, becomes with , are two integration constants, and is the normalized phase velocity. If has real roots, (2.54) has the stationary solution where is Jacobi’s elliptic function of modulus [39], is the amplitude, is a constant of integration, andare the complete elliptic integrals [40].
Parameter represents measure characterizing the degree of nonlinearity in the system. Therefore, the solution (2.57) contains (as subsequences for ) one-dimensional harmonic waves, while for, one-dimensional wave packet. These two subsequences define the nonquasiautonomous regime of the drug release process [22, 36, 37], that is, the system should receive external energy in order to develop. For , the solution (2.57) becomes one-dimensional soliton, while for , one-dimensional soliton packet will be generated. The last two imply a quasiautonomous regime (self-evolving and independent [22]) for drug particle release process [22, 36, 37].
The three-dimensional plot of solution (2.57) shows one-dimensional cnoidal oscillation modes of the concentration field, generated by similar trajectories of the drug particles (see Figure 3). We mention that cnoidal oscillations are nonlinear ones, being described by the elliptic function cn, hence the name (cnoidal).
It is known that in nonlinear dynamics, cnoidal oscillation modes are associated with nonlinear lattice of oscillators (the Toda lattice [41]). Consequently, large-time-scale drug particle ensembles can be compared to a lattice of nonlinear oscillators which facilitates drug release process.
3. Experimental Results
Most of the experimental data in the literature reveal that, on average, drug release from polymeric matrices takes place according to a power law in the first 60% of the release curve and/or to exponential Weibull function on the entire drug release curve, reaching an average constant balanced value. The majority of these experimental results are carried out on relatively short time intervals, dissolution and diffusion being the dominant systems. The system exhibits a “burst effect” due to highly concentrated gradient. The phenomenon is followed by linear evolution on a constant value that corresponds to the balanced state with diluted gradient.
Nevertheless, some experimental results with long enough time intervals allow complete evolution of the process (polymer degradation stage is included here) and show unusually strong fluctuating behaviour.
Experimental data of drug release, at short and long time scales, for polymeric microparticles (as polymeric matrices) are presented below.
3.1. Materials
The following materials were used: low-molecular-weight chitosan-CS, deacetylation degree 75–80% (Aldrich), type B gelatin-GEL (Aldrich), glutaraldehyde-GA (Aldrich)-25% aqueous solution, sodium tripolyphosphate-TPP (Sigma), Levofloxacin-LEV (Sigma), Tween 80 (Aldrich), and Span 80 (Aldrich).
3.2. Preparation of Microparticles by Ionic Gelation and Covalent Cross-Linking in O/W/O Emulsion
Microparticles were prepared using an original double cross-linking method of a CS-GEL mixture. Different weight ratios CS/GEL (in terms of amino groups of both polymers) were dissolved in acetic acid solution 2%, and then Tween 80 was added to make a 2% (w/w) surfactant in the solution. The mixture was magnetically stirred until the surfactant was completely dissolved. Two different solutions of 2% Span 80 in toluene were prepared according to O1/W (v/v) = 1/4 and O2/(O1 + W) = 4/1. The organic phase O1 was dripped within the aqueous polymer phase, W under homogenization with an Ultraturax device at 9000 rpm. The primary emulsion was transformed into a double one through dripping in the second organic phase O2, according to the same hydrodynamic regime. The emulsion was then gelled by slowly adding a TPP solution at a rate of 2 mL/min with continuous stir for extra 10 min.
The suspension was then transferred to a round-bottom flask and mechanically stirred at 500 rpm. A certain amount of a saturated solution of GA in toluene was consecutively added and stirred for 60 min. The particles were separated by centrifugation (6000 rpm) and repeatedly washed with acetone and water in order to eliminate residual compounds. After hexane wash, the particles were dried at room temperature.
3.3. Preparation Parameters
A two-step solidification method was used. The first step, which has critical influence over the subsequent particle shape and size, included ionic cross-linking with TPP effect through phosphate bridges among amino functionalities in both types of polymeric chains. The GA covalent cross-linking (also taking place in NH2 groups) was performed with the purpose of stabilizing gel capsules. Our study analyzes the influence of the following cross-linking reaction parameters on the levofloxacin release kinetics:(i)concentration of the ionic cross-linker,(ii)ratio among amino functionalities of the two polymers and the ionic crosslinker,(iii)polymer composition of the polymer mixture.
Table 1 shows the variable parameters in the preparation protocol that have been grouped according to the variable parameters.
3.4. Levofloxacin Release Kinetics
3.4.1. Levofloxacin Release Kinetics at Small-Time Scales
If the experimental time scale is of minutes order, the evolution of the released drug concentration will be described by Peppas law. In this case, the correlated factor ranges between 0.8413 and 0.9983. Experimental and Peppas curves can be observed in Figure 4 (the Peppas parameters and the correlation coefficient for each sample are given in Table 2). The plots group according to variable preparation parameters (for a better observation of the first points, time scale is 500 min, although the fitting was made on the points up to 1440 min (one day)). Relative errors range between 1% and 5%, with no important influence on release kinetic evolution.
(a)
(b)
(c)
Previous works have shown the form dependence [38, 39] between the value of parameter in Peppas equation (considered as short-time approximation of Weibull function) and the fractal dimension of the drug particle during the release process ()
Thus, according to experimental data, the following values were obtained in Table 2.
One first observation refers to the proportional dependence among experimental variable, on one hand, and Peppas parameters, on the other, in the particular case of the third sample group, that is, increases with the chitosan/gelatin ratio. This proves to be experimentally useful if we want to obtain a Fickian diffusion. At the same time, the concentration of the released drug proves to be very low. This could be explained by drug crystallization inside the microparticle and the dependence of its release (dissolution followed by diffusion) on polymer degradation.
In our opinion, the value of the fractal dimension is important as long as its values are unusually high and indicate that either fractal dimension must be considered as function of structure “classes,” or drug release processes (implicitly drug particle trajectories) have high degrees of complexity and nonlinearities, implying many freedom degrees in the phase space [42].
This analysis (small concentration of the released drug and high fractal dimensions) made us continue the experiment until the system reached a stationary state.
3.4.2. Levofloxacin Release Kinetics at Large-Time Scales
The experiments at large-time scales (of days order) revealed unusual behavior characterized by large variations. The release kinetics of levofloxacin is plotted in Figure 5. The relative errors range between 1% and 5%, and, for better visualization, the error bars are plotted in Figure 6.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Experiments have been performed for 28 days, the concentration of the released drug being measured daily, at the same hour. The general characteristic of the above kinetics refers to strong variations of concentration in time, approximately at the same moment.
In the following section, we will explain the evolution of these systems through the theoretical model (developed in Section 2) based on fractal approximation of motion.
3.5. The Correspondence between Theoretical Model and Experimental Results
In what follows, we identify the field from relation (2.57) with normalized concentration field of the released drug from microparticles.
For best correlation between experimental data and the theoretical model (for each sample), we used a planar intersection of the graph in Figure 3, where the two variables are and . With these variables, (2.57) becomes
Thus, in order to find the one-dimensional equation for a planar intersection, perpendicular to plane , we used (linear function equation), where and are two parameters. This equation is transformed into a parametric equation by means of the following substitutions:in (3.2).
Afterwards, we obtain one-dimensional function
The highest value of the correlation coefficient (for two vectors: one obtained from this very function, the other from experimental data) for different values of and (in the particular experimental case) will represent the best approximation of experimental data with the theoretical model.
Our goal was to find the right correlation coefficient which should be higher than , in order to demonstrate the relevance of the model we had in view. Figure 6 shows experimental and theoretical curves that were obtained through our method, where represents the correlation coefficient and a normalized variable which is simultaneously dependent on normalized time and on nonlinear degree of the system ( parameter). Geometrically, represents the congruent angle formed by the time axis and the vertical intersection plane.
Parameters and of the planar intersections for the above theoretical curves are shown in Table 3.
We must mention that for each sample the fitting process was an independent one. The corresponding intersection plane that offers best correlation factors had to be identified by each sample in turn.
A first observation refers to the correlation among plane and variable parameters (within experimental protocol) differ from Peppas small-time-scale fitting.
We consider that this could be a starting point in establishing dependence among experimental parameters involved in the protocol. The purpose of this analysis is to obtain polymer matrices together with characteristics of release kinetics, taking into account that until now, this type of dependence had to pass through intermediary stages of the physical and chemical characterization of polymer matrices.
The few experimental data could not sustain a general conclusion on the existing dependence among plane and experimental parameters, but this will be the purpose of a next paper.
4. Conclusions
If the particle moves on fractal curves, a new model for drug release mechanism from polymer matrix (namely, polymeric particles) is obtained. This model offers new alternatives for the theoretical study of drug release process (on large time scale) in the presence of all phenomena and considering a highly complex and implicitly nonlinear system. Consequently, the concentration field has cnoidal oscillation modes, generated by similar trajectories of drug particles. This means that the drug particle ensemble (at time large scale) works in a network of nonlinear oscillators, with oscillations around release boundary. Moreover, the normalized concentration field simultaneously depends on normalized time nonlinear system (through parameter).
The fitting procedure among experimental and theoretical curves revealed the existing correlation of some characteristics of the release kinetics (the parameters of the intersection plane) with variable experimental parameters.
At the same time, we consider that this could be a starting point in establishing dependence among experimental parameters, taking into account that until now, this type of dependence had to pass through intermediary stages of physical and chemical characteristics of polymer matrices.
Acknowledgments
This paper was supported by the project PERFORM-ERA “Postdoctoral Performance for Integration in the European Research Area” (ID-57649), financed by the European Social Fund and the Romanian Government, and by the project POSDRU/88/1.5/S/47646 of The European Social Fund.