Nonlinearities in Drug Release Process from Polymeric Microparticles: Long-Time-Scale Behaviour

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 𝑥_𝑎=𝑥₀,𝑥₁=𝑥₀+𝜀,…,𝑥_𝑘=𝑥₀+𝑘𝜀,𝑥_𝑛=𝑥₀+𝑛𝜀=𝑥_𝑏.(2.1)

Let us note by 𝐹(𝑥,𝜀) the broken line that connects the points 𝐹𝑥₀,…,𝐹𝑥_𝑘,…,𝐹𝑥_𝑛.(2.2)

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 𝑑𝜀𝜀=𝑑𝜀𝜀=𝑑𝜌,(2.3) a situation in which we can consider the infinitesimal-scale transformation as being 𝜀^=𝜀+𝑑𝜀=𝜀+𝜀𝑑𝜌.(2.4)

Figure 1 

Dilatation-scale invariance.

Such transformation in the case of function 𝐹(𝑥,𝜀) leads to 𝐹𝑥,𝜀^=𝐹(𝑥,𝜀+𝜀𝑑𝜌),(2.5) respectively, if we limit ourselves to a first-order approximation: 𝐹𝑥,𝜀^=𝐹(𝑥,𝜀)+𝜕𝐹(𝑥,𝜀)𝜕𝜀(𝜀′−𝜀)=𝐹(𝑥,𝜀)+𝜕𝐹𝜕(𝑥𝜀,𝜀)𝜀𝑑𝜌.(2.6)

Moreover, let us notice that for an arbitrary but fixed 𝜀₀, we obtain 𝜕ln𝜀/𝜀₀𝜕𝜀=𝜕ln𝜀−ln𝜀⁰𝜕𝜀=1𝜀,(2.7) a situation in which (2.6) can be written as follows: 𝐹𝑥,𝜀^=𝐹(𝑥,𝜀)+𝜕𝐹(𝑥,𝜀)𝜕ln𝜀/𝜀₀𝑑𝜌=1+𝜕ln𝜕𝜀/𝜀₀𝑑𝜌𝐹(𝑥,𝜀).(2.8)

Therefore, we can introduce the dilatation operator 𝐷=𝜕𝜕ln𝜀/𝜀₀.(2.9)

At the same time, relation (2.9) shows that the intrinsic variable of resolution is not 𝜀, but ln(𝜀/𝜀₀).

The fractal function is explicitly dependant on the resolution (𝜀/𝜀₀); therefore, we have to solve the differential equation 𝑑𝐹𝑑ln𝜀/𝜀₀=𝑃(𝐹),(2.10) where 𝑃(𝐹) is now an unknown function. The simplest explicit suggested form for 𝑃(𝐹) is linear dependence [29] 𝑃(𝐹)=𝐴+𝐵𝐹,𝐴,𝐵=const.,(2.11) in which case the differential equation (2.10) takes the form 𝑑𝐹𝑑ln𝜀/𝜀₀=𝐴+𝐵𝐹.(2.12)

Hence, by integration and substituting−𝐵𝐴=−𝜏,(2.13a)𝐵=𝐹⁰,(2.13b)we obtain 𝐹𝜀𝜀₀=𝐹⁰1+𝜀𝜀⁰^𝜏.(2.14)

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 𝜀/𝜀₀) 𝐹𝑝,𝜀𝜀₀=𝐹⁰(𝑝)1+𝜉(𝑝)𝜀⁰𝜀^𝜏(𝑝),(2.15) where 𝜉(𝑝) is a function depending on parameterization of the fractal curve.

The following particular cases are to be considered.

(ii1) In asymptotic small-scale regime𝜀≪𝜀₀, 𝜏 is constant (with no scale dependence) and power-law dependence on resolution is obtained:𝐹𝑝,𝜀𝜀₀=𝑇(𝑝)𝜀𝜀⁰^𝜏𝑇(𝑝)=𝐹₀,(2.16a)(𝑝)𝑄(𝑝).(2.16b)

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: 𝑀_𝑡𝑀_∞=𝑘_𝐻⋅𝑡^1/2,(2.17)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: 𝑀_𝑡𝑀_∞=𝑘⋅𝑡^𝑛,(2.18)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, 𝑛=0.5 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 𝐷_𝑇=1, since we deal with length, but this can be easily generalized to surfaces (𝐷_𝑇=2) and volumes (𝐷_𝑇=3). 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 𝑋(𝑝,𝛿𝑋)=𝑋₀(𝑝)⋅𝜆𝛿𝑋^𝐷𝐹−1,(2.19) where 𝜆 is a scale characteristic length, and the exponent is identified with 𝜏=𝐷_𝐹−1.

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 𝛿𝑋=𝑎𝛿𝑡⋅𝜆𝛿𝑋^𝐷𝐹−1,(2.20) so that the following fundamental relation among space and time elements on a fractal curve or function is obtained: 𝛿𝑋^𝐷𝐹∝𝛿𝑡.(2.21) In other words, they are differential elements of different orders.

(ii2) In the asymptotic big-scale regime 𝜀≫𝜀₀, 𝜏 is constant (with no scale dependence), and, in terms of resolution, one obtains an independent law 𝐹𝑝,𝜀𝜀₀→𝐹⁰(𝑝).(2.22)

Particularly, if 𝐹(𝑝,𝜀/𝜀₀) are the coordinates in given space, we can write 𝑋𝑝,𝜀𝜀₀=𝑥(𝑝)1+𝜉(𝑝)𝜀𝜀⁰^𝜏.(2.23)

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 𝑑𝑋=𝑑𝑥+𝑑𝜉,(2.24) 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:𝑑₊𝑄𝑑𝑡=lim^Δ𝑡→0₊𝑄(𝑡+Δ𝑡)−𝑄(𝑡)𝑑Δ𝑡,(2.25a)₋𝑄𝑑𝑡=lim^Δ𝑡→0₋𝑄(𝑡)−𝑄(𝑡−Δ𝑡)Δ𝑡.(2.25b)

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 𝑑𝑡→0. As a consequence, this leads to the two derivatives of the fractal field 𝑄 as explicit functions of the two variables 𝑡 and 𝑑𝑡,𝑑₊𝑄𝑑𝑡=lim^Δ𝑡→0₊𝑄(𝑡+Δ𝑡,Δ𝑡)−𝑄(𝑡,Δ𝑡)𝑑Δ𝑡,(2.26a)₋𝑄𝑑𝑡=lim^Δ𝑡→0₋𝑄(𝑡,Δ𝑡)−𝑄(𝑡−Δ𝑡,Δ𝑡)Δ𝑡.(2.26b)

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.

Figure 2 

The continuous curves which are not fractals but have certain points where they are not differentiable.

If we consider all the lines (segments) starting from 𝑃, we denote the average of these vectors by 𝑑𝑥^𝑖±, that is, 𝑑𝑑𝑋𝑋^𝑖+_𝑖−=𝑑𝑥^𝑖+,𝑖=1,2,(2.27a)=𝑑𝑥^𝑖−,𝑖=1,2.(2.27b)

Since, according to (2.24), we can write𝑑𝑑𝑋𝑋^𝑖+_𝑖−=𝑑𝑥^𝑖++𝑑𝜉^𝑖+,(2.28a)=𝑑𝑥^𝑖−+𝑑𝜉^𝑖−,(2.28b)

it results that𝑑𝑑𝜉𝜉^𝑖+_𝑖−=0,(2.29a)=0.(2.29b)

(vi) The differential fractal part satisfies, according to (2.21), the fractal equation𝑑𝑑₊𝜉^𝑖=𝜆^𝑖+(𝑑𝑡)^1/𝐷𝐹,(2.30a)₋𝜉^𝑖=𝜆^𝑖−(𝑑𝑡)^1/𝐷𝐹,(2.30b)

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 𝑑𝑑𝑡=12𝑑⁺+𝑑₋𝑑𝑡−2𝑖𝑑⁺−𝑑₋𝑑𝑡.(2.31)

Applying this operator to the “position vector,” a complex speed yields 𝐕=𝑑𝐗𝑑𝑡=12𝑑⁺𝐗+𝑑₋𝐗𝑑𝑡−2𝑖𝑑⁺𝐗−𝑑₋𝐗𝑑𝑡=𝐕⁺+𝐕₋2−𝑖𝐕⁺−𝐕₋2=𝐕−𝑖𝐔,(2.32) with𝐕=𝐕⁺+𝐕₋𝐔=𝐕₊2,(2.33a)−𝐕₋2.(2.33b)

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 𝑋^𝑖(𝑖=1,3) 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 order𝑑₊𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+∇𝑄⋅𝑑⁺𝐗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₊𝑋^𝑖𝑑₊𝑋^𝑗+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑𝑑₊𝑋^𝑖𝑑₊𝑋^𝑗𝑑₊𝑋^𝑘,(2.34a)₋𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+∇𝑄⋅𝑑⁻𝐗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₋𝑋^𝑖𝑑₋𝑋^𝑗+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑₋𝑋^𝑖𝑑₋𝑋^𝑗𝑑₋𝑋^𝑘,(2.34b)where 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 𝑑_±𝑄=𝜕𝜕𝑄𝑡𝑑𝑡+∇𝑄⋅𝑑^±𝑋+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑_±𝑋^𝑖𝑑_±𝑋^𝑗+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑_±𝑋^𝑖𝑑_±𝑋^𝑗𝑑_±𝑋^𝑘⟨𝑑₋,(2.35)𝑄⟩=𝜕𝑄𝜕𝑡𝑑𝑡+⟨∇𝑄⋅𝑑⁻𝑋⟩+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₋𝑋^𝑖𝑑₋𝑋^𝑗+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑₋𝑋^𝑖𝑑₋𝑋^𝑗𝑑₋𝑋^𝑘.(2.36)

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𝑑₊𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+∇𝑄𝑑⁺𝐗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₊𝐗^𝑖𝑑₊𝐗^𝑗+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑𝑑₊𝐗^𝑖𝑑₊𝐗^𝑗𝑑₊𝐗^𝑘,(2.37a)₋𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+∇𝑄⟨𝑑⁻𝐗⟩+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₋𝐗^𝑖𝑑₋𝐗^𝑗+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑₋𝐗^𝑖𝑑₋𝐗^𝑗𝑑₋𝐗^𝑘,(2.37b)

or more, using (2.28a) and (2.28b) with characteristics (2.29a) and (2.29b),𝑑₊𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+∇𝑄⋅𝑑⁺𝐗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+1𝑑₊𝐱^𝑖𝑑₊𝐱^𝑗+𝑑₊𝜉^𝑖𝑑₊𝜉^𝑗6𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑𝑑₊𝐱^𝑖𝑑₊𝐱^𝑗𝑑₊𝐱^𝑘+𝑑₊𝜉^𝑖𝑑₊𝜉^𝑗𝑑₊𝜉^𝑘,(2.38a)₋𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+∇𝑄⋅𝑑⁻𝐗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+1𝑑₋𝐱^𝑖𝑑₋𝐱^𝑗+𝑑₋𝜉^𝑖𝑑₋𝜉^𝑗6𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑₋𝐱^𝑖𝑑₋𝐱^𝑗𝑑₋𝐱^𝑘+𝑑₋𝜉^𝑖𝑑₋𝜉^𝑗𝑑₋𝜉^𝑘.(2.38b)

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𝑑𝑑₊𝜉^𝑖𝑑₊𝜉^𝑗=𝜆^𝑖+𝜆^𝑗+(𝑑𝑡)^{(2/𝐷𝐹)−1}𝑑𝑡,(2.39a)₋𝜉^𝑖𝑑₋𝜉^𝑗=𝜆^𝑖−𝜆^𝑗−(𝑑𝑡)^{(2/𝐷𝐹)−1}𝑑𝑡.(2.39b)

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𝑑𝑑₊𝜉^𝑖𝑑₊𝜉^𝑗𝑑₊𝜉^𝑘=𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+(𝑑𝑡)^{(3/𝐷𝐹)−1}𝑑𝑡,(2.40a)₋𝜉^𝑖𝑑₋𝜉^𝑗𝑑₋𝜉^𝑘=𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−(𝑑𝑡)^{(3/𝐷𝐹)−1}𝑑𝑡.(2.40b)

Then, (2.38a) and (2.38b) may be written as follows:𝑑₊𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+𝑑⁺𝐱⋅∇𝑄+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₊𝐱^𝑖𝑑₊𝐱^𝑗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+1𝜆^𝑖+𝜆^𝑗+(𝑑𝑡)^{(2/𝐷𝐹)−1}𝑑𝑡6𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑₊𝐱^𝑖𝑑₊𝐱^𝑗𝑑₊𝐱^𝑘+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+(𝑑𝑡)^{(3/𝐷𝐹)−1}𝑑𝑡,(2.41a)₋𝑄=𝜕𝑄𝜕𝑡𝑑𝑡+𝑑⁻𝐱⋅∇𝑄+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝑑₋𝐱^𝑖𝑑₋𝐱^𝑗+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+1𝜆^𝑖−𝜆^𝑗−(𝑑𝑡)^{(2/𝐷𝐹)−1}𝑑𝑡6𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑₋𝐱^𝑖𝑑₋𝐱^𝑗𝑑₋𝐱^𝑘+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−(𝑑𝑡)^{(3/𝐷𝐹)−1}𝑑𝑡.(2.41b)

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𝑑₊𝑄𝑑𝑡=𝜕𝜕𝑄𝑡+𝐕⁺⋅∇𝑄+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜆^𝑖+𝜆^𝑗+(𝑑𝑡)^{(2/𝐷𝐹)−1}+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝑑𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+(𝑑𝑡)^{(3/𝐷𝐹)−1},(2.42a)₋𝑄𝑑𝑡=𝜕𝜕𝑄𝑡+𝐕⁻⋅∇𝑄+12𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜆^𝑖−𝜆^𝑗−(𝑑𝑡)^{(2/𝐷𝐹)−1}+16𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−(𝑑𝑡)^{(3/𝐷𝐹)−1}.(2.42b)

These relations also allow us to define the operator𝑑₊𝑑𝑡=𝜕𝜕𝑡+𝐕⁺⋅∇+12𝜕𝑋𝜕^𝑖𝜕²𝑋^𝑗𝜆^𝑖+𝜆^𝑗+(𝑑𝑡)^{(2/𝐷𝐹)−1}+16𝜕𝑋^𝑖𝜕𝜕𝑋^3𝑗𝜕𝑋^𝑘𝑑𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+(𝑑𝑡)^{(3/𝐷𝐹)−1},(2.43a)₋𝑑𝑡=𝜕𝜕𝑡+𝐕⁻⋅∇+12𝜕𝑋𝜕^𝑖𝜕²𝑋^𝑗𝜆^𝑖−𝜆^𝑗−(𝑑𝑡)^{(2/𝐷𝐹)−1}+16𝜕𝑋^𝑖𝜕𝜕𝑋^3𝑗𝜕𝑋^𝑘𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−(𝑑𝑡)^{(3/𝐷𝐹)−1}.(2.43b)

Under these circumstances, let us calculate (𝜕̂𝑄/𝜕𝑡). Taking into account (2.43a), (2.43b), (2.31), and (2.32), we will obtain 𝜕^∧𝑄𝜕𝑡=12𝑑⁺𝑄𝑑𝑡+𝑑⁻𝑄𝑑𝑡−𝑖𝑑⁺𝑄𝑑𝑡−𝑑⁻𝑄=1𝑑𝑡2𝜕𝜕𝑄𝑡+12𝐕⁺⋅∇𝑄+𝜆^𝑖+𝜆^𝑗+14(𝑑𝑡)^{(2/𝐷𝐹)−1}𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+112(𝑑𝑡)^{(3/𝐷𝐹)−1}𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘+12𝜕𝜕𝑄𝑡+12𝐕⁻⋅∇𝑄+𝜆^𝑖−𝜆^𝑗−14(𝑑𝑡)^{(2/𝐷𝐹)−1}𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−112(𝑑𝑡)^{(3/𝐷𝐹)−1}𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘−𝑖2𝜕𝜕𝑄𝑡−2𝑖𝐕⁺⋅∇𝑄−𝜆^𝑖+𝜆^𝑗+𝑖2(𝑑𝑡)^{(2/𝐷𝐹)−1}𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗−𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+𝑖12(𝑑𝑡)^{(3/𝐷𝐹)−1}𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘+𝑖2𝜕𝜕𝑄𝑡+2𝑖𝐕⁻⋅∇𝑄+𝜆^𝑖−𝜆^𝑗−𝑖2(𝑑𝑡)^{(2/𝐷𝐹)−1}𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−𝑖12(𝑑𝑡)^{(3/𝐷𝐹)−1}𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘=𝜕𝑄𝜕𝑡+𝐕⁺+𝐕₋2−𝑖𝐕⁺−𝐕₋+(𝑑𝑡)_{(2/𝐷𝐹)−1}2⋅∇𝑄4𝜆^𝑖+𝜆^𝑗++𝜆^𝑖−𝜆^𝑗−−𝑖𝜆^𝑖+𝜆^𝑗+−𝜆^𝑖−𝜆^𝑗−𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+(𝑑𝑡)^{(3/𝐷𝐹)−1}12𝜆^𝑖+𝜆^𝑗+𝜆^𝑘++𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−−𝑖𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+−𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘=𝜕𝑄𝜕𝑡+^∧𝐕⋅∇𝑄+(𝑑𝑡)^{(2/𝐷𝐹)−1}4𝜆^𝑖+𝜆^𝑗++𝜆^𝑖−𝜆^𝑗−−𝑖𝜆^𝑖+𝜆^𝑗+−𝜆^𝑖−𝜆^𝑗−𝜕²𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗+(𝑑𝑡)^{(3/𝐷𝐹)−1}12𝜆^𝑖+𝜆^𝑗+𝜆^𝑘++𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−−𝑖𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+−𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−𝜕³𝑄𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘.(2.44)

This relation also allows us to define the fractal operator 𝜕̂𝜕𝑡=𝜕𝜕𝑡+𝐕⋅∇+(𝑑𝑡)4^{(2/𝐷𝐹)−1}𝜆^𝑖+𝜆^𝑗++𝜆^𝑖−𝜆^𝑗−−𝑖𝜆^𝑖+𝜆^𝑗+−𝜆^𝑖−𝜆^𝑗−𝜕²𝜕𝑋^𝑖𝜕𝑋^𝑗+(𝑑𝑡)^{(3/𝐷𝐹)−1}12𝜆^𝑖+𝜆^𝑗+𝜆^𝑘++𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−−𝑖𝜆^𝑖+𝜆^𝑗+𝜆^𝑘+−𝜆^𝑖−𝜆^𝑗−𝜆^𝑘−𝜕³𝜕𝑋^𝑖𝜕𝑋^𝑗𝜕𝑋^𝑘.(2.45)

Particularly, by choosing𝜆𝜆^𝑖+𝜆^𝑗+=−𝜆^𝑖−𝜆^𝑗−=2𝒟𝛿^𝑖𝑗,(2.46a)_𝑖+𝜆^𝑗+𝜆^𝑘+=−𝜆^𝑖−𝜆^𝑗−𝜆^𝑘+=2√2𝒟^3/2𝛿^𝑖𝑗𝑘,(2.46b)the fractal operator (2.45) takes the usual form 𝜕̂𝜕𝑡=𝜕𝜕𝑡+𝐕⋅∇−𝑖𝒟(𝑑𝑡)^{(2/𝐷𝐹)−1}Δ+√23𝒟^3/2(𝑑𝑡)^{(3/𝐷𝐹)−1}∇³.(2.47)

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 𝜕̂𝑄𝜕𝑡=𝜕𝜕𝑄𝑡+𝑉⋅∇𝑄−𝑖𝒟(𝑑𝑡)^{(2/𝐷𝐹)−1}Δ𝑄+√23𝒟^3/2(𝑑𝑡)^{(3/𝐷𝐹)−1}∇³𝑄=0.(2.48)

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 𝜕̂𝑄𝑑𝑡=𝜕𝜕𝑄𝑡+𝐕⋅∇𝑄+√23𝐷^3/2(𝑑𝑡)^{(3/𝐷𝐹𝐷)−1}∇³𝑄=0.(2.49)

If we separate the real and imaginary parts from (2.49), we will obtain𝜕𝑄𝜕𝑡+𝐕⋅∇𝑄+√23𝐷^3/2−𝐔⋅∇𝑄(=𝑑𝑡0).^{(3/𝐷𝐹)−1}∇³𝑄=0,((22..5500ba))

By adding them, the fractal diffusion equation is 𝜕𝑄𝜕𝑡+(𝐕−𝐔)⋅∇𝑄+√23𝐷^3/2(𝑑𝑡)^{(3/𝐷𝐹)−1}∇³𝑄=0.(2.51)

From (2.50b), we see that, at fractal scale, there will be no 𝑄 field gradient.

Assuming that |𝑉−𝑈|=𝜎⋅𝑄 with 𝜎=constant (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𝜏=𝜔𝑡,(2.52a)Φ𝜉==𝑘𝑄𝑥,(2.52b)𝑄₀,(2.52c)

and normalizing conditions 𝜎𝑄₀𝑘6𝜔=√23𝐷^3/2(𝑑𝑡)^{(3/𝐷𝐹)−1}𝑘³𝜔=1(2.53) take the form 𝜕_𝜏𝜙+6𝜙𝜕_𝜉𝜙+𝜕_𝜉𝜉𝜉𝜙=0.(2.54)

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.55a)𝜉−𝑢𝜏,(2.55b)

(2.54), by double integration, becomes 12𝑤^2=𝐹(𝑤)=−𝑤³−𝑢2𝑤²−𝑔𝑤−ℎ,(2.56) with 𝑔, ℎ are two integration constants, and 𝑢 is the normalized phase velocity. If 𝐹(𝑤) has real roots, (2.54) has the stationary solution 𝜙𝜉,𝜏,𝑠=2𝑎𝐸(𝑠)𝐾(𝑠)−1+2𝑎⋅cn²√𝑎𝑠𝜉−𝑢2𝜏+𝜉₀;𝑠,(2.57) where cn is Jacobi’s elliptic function of 𝑠 modulus [39], 𝑎 is the amplitude, 𝜉₀ is a constant of integration, and𝐾𝐸((𝑠𝑠))==^0𝜋/2_0𝜋/21−𝑠²sin²𝜑^−1/2𝑑𝜑,(2.58a)1−𝑠²sin²𝜑^1/2𝑑𝜑,(2.58b)are 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 𝑠=0) one-dimensional harmonic waves, while for, 𝑠→0 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 𝑠=1, the solution (2.57) becomes one-dimensional soliton, while for 𝑠→1, 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).

Figure 3 

One-dimensional cnoidal oscillation modes of the field Φ.

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 NH₂ 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.

Table 1 

The variable parameters in the preparation protocol.

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.

Table 2 

The variable parameters of the preparation protocol.

(a)

(b)

(c)

Figure 4 

Levofloxacin release kinetics (experimental and Peppas fitting), at time scale of hours order, at different concentration of TPP solution (a), NH2: TPP mol ratios (b), and NH2 mols (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 (𝐷_𝑓) 𝑛=2𝐷_𝑓.(3.1)

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)

Figure 5 

Experimental release kinetics of levofloxacin, at time scale of days order, at different concentration of TPP solution (a), NH2: TPP mol ratios (b), and NH2 mols (c).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 6 

The best correlations among experimental and theoretical curves (blue line—experimental curve, red line—theoretical curve).

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 𝑦=(𝜉−𝜏𝑢)/2 and 𝑥=𝑠. With these variables, (2.57) becomes𝜙₁(𝑥,𝑦)=2𝑎𝐸(𝑥)𝐾(𝑥)−1+2𝑎⋅cn²√𝑎𝑥𝑦+𝜉⁰;𝑥.(3.2)

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:𝑥=𝑙𝑦=𝑛+√𝑚𝑚²+𝑙1,(3.3a)√𝑚²+1,(3.3b)in (3.2).

Afterwards, we obtain one-dimensional function𝜙₂(𝑡,𝑚,𝑛)=𝜙₁√𝑚²𝑡+1,𝑛+𝑚√𝑚²𝑡+1.(3.4)

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 0.6-0.7, 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.

Table 3 

Parameters of the planar intersections.

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.