A theoretical mathematical model for assessing diclofenac release from chitosan-based formulations

Abstract The paper reports a new mathematical model for understanding the mechanism delivery from drug release systems. To do this, two drug release systems based on chitosan and diclofenac sodium salt as a drug model, were prepared by in situ hydrogelation in the presence of salicylaldehyde. The morphology of the systems was analyzed by scanning electron microscopy and polarized light microscopy and the drug release was in vitro investigated into a medium mimicking the in vivo environment. The drug release mechanism was firstly assessed by fitting the in vitro release data on five traditional mathematical model. In the context of pharmacokinetics behavioral analysis, a new mathematical procedure for describing drug release dynamics in polymer-drug complex systems was proposed. Assuming that the dynamics of polymer-drug system’s structural units take place on continuous and nondifferentiable curves (multifractal curves), it was showed that in a one-dimensional hydrodynamic formalism of multifractal variables the drug release mechanism is given through synchronous dynamics at a differentiable and non-differentiable scale resolutions.


Introduction
The rapid technological development contributed on a hand to the improvement of the human life and on the other hand to the increase of the mortal illnesses impact, such as tumors, diabetes and heart attack, and so on. The dramatic side effects of the systemic drug administration directed the research to a new domain, that of the controlled drug delivery systems, which become an important topic of our days (Tiwari et al., 2012). Basically, these systems consist in the dispersion of the drug into a matrix, which has the ability to deliver it in a controlled, prolonged manner to a targeted site. From this perspective, the development of the domain implies the investigation of new systems and of the mechanisms which governs the drug release. There are a lot of drug delivery systems explored in the recent times, but there are still certain trials that must be addressed for a successful delivery of the drugs to target site (Patra et al., 2018;Fisher et al., 2010). To be used as drug delivery systems in medical field, the corresponding matrix needs to present specific properties such as biocompatibility, biodegradability and ability to develop strong forces with the drug in order to prolong its release (Peppas & Brannon-Peppas, 2001;Craciun et al., 2019). Different materials, derived from natural or synthetic resources have been used as matrix for drug delivery systems over time (Griesser & Kambouris, 2001). By far the hydrogels, especially those obtained from natural resources, appeared to meet the requisite requirements for drugs matrix, due to the fact that besides biocompatibility, biodegradability, and lack of toxicity of the degradation byproducts, they also have the advantage to mimic tissues by a proper moisture degree. Chitosan based hydrogels successfully fulfills all these required properties, being preferred for applications in medicine, but also in agriculture, water waste treatment, hygiene, food industry, and so on (Ahmadi et al., 2015;Desbri eres & Guibal, 2018;Iftime et al., 2019;Abdel-Aleem & El-Aidie, 2018). The properties of the chitosan hydrogels are closely related with their preparation pathway. Those prepared by physical crosslinking are able to swell and load different drugs, but most of them have a too fast degradation rate over time (Vitaliy & Khutoryanskiy, 2015). On the other hand, those prepared by chemical cross-linking have a better stability but their degradation rate is often too low for in vivo applications (Peppas et al., 2000). Besides, an important issue is the chemical crosslinker, which is usually a dialdehyde with a certain degree of toxicity limiting the hydrogel bio-application (Omidian & Park, 2010). Taking account of this state of the art, our group developed a new strategy for chitosan hydrogelation with natural monoaldehydes, which are of natural origin and have no toxicity (Marin et al., 2015;Ailincai et al., 2016;Iftime et al, 2017;Marin et al., 2017;Olaru et al., 2018;Bejan et al., 2018;Iftime & Marin, 2018;Xiong et al., 2019). The hydrogels prepared by this route proved excellent ability to accommodate biologic agents, such as drugs or fertilizers (Craciun et al., 2019;Iftime et al., 2019;Ailincai et al., 2018;Iftime et al., 2020). Among these new hydrogels, that using salicylaldehyde as crosslinker lead to chitosan hydrogels which demonstrated appropriate properties for the design of drug delivery systems, such as biocompatibility, biodegradability, swelling ability, thixotropy and self-healing (Iftime et al, 2017). They proved capacity to encapsulate a model drug, diclofenac sodium salt (DCF), and further to slowly release it along 10 days, based on the strong physical forces developed between matrix and drug (Iftime et al., 2020). These excellent findings, encouraged us to further investigate these formulations in order to better understand the release mechanism, to create a proper structure-drug release correlation, and further to improve the design to reach the thorny goal of a real world application. Since the drug release is a complex phenomenon which depends by many factors, in this paper for a better understanding of the DCF release mechanism from the chitosan based hydrogels it was proposed to model the drug release kinetics by a multifractal theoretical model built by means of empiric laws and logistics type. Traditional mathematic models, such as Zero order, First order, Korsmeyer-Peppas, Higuchi model and Hixson-Crowell model were used in order to validate the multifractal model developed by us.
The above mentioned models are based on the hypothesis of homogeneity in its various forms (homogenous kinetic space, law of mass etc.) which has become almost dogmatic in classical Pharmacokinetics (PK). The functionality of such a hypothesis allowed the development of a class of differentiable models in the description of dynamics of biological systems (i.e., 'compartmental' analysis) and mainly, of drug release dynamics in such systems. However, biological systems are nowadays understood as inherently nondifferential (fractal). Specifically, the microenvironments where any drug molecules with membrane interface, metabolic enzymes or pharmacological receptors are unanimously recognized as unstirred, spacerestricted, heterogeneous and geometrically fractal. It is thus necessary to define a new class of models, this time nondifferentiable, in describing biological system dynamics and particularly drug release dynamics in such systems. The Fractal Pharmacokinetics implies the use of fractional calculus, expanding on the notion of dimension etc. As such, it is possible in the context of 'compartmental analysis' (Pereira, 2010) to describe diffusion in dense objects (Lemehaute & Crepy, 1983), dynamics in polymeric networks (Barkai & Klafter, 1998), diffusion in porous and fractal media (O'Shaughnessy & Procaccia, 1985), kinetics in viscoelastic media (Mainardi, 1994) etc. More recently, 'compartmental analysis' through PK allowed the modeling of processes such as drug dissolution (Kosmidis et al., 2003), absorption (Higaki et al, 2001), distribution (Karalis et al., 2003), whole disposition (Weiss, 1999), kinetics with biomolecular reactions (Kotulskil & Weron, 1996) etc. In this work a new method for describing drug release dynamics in complex systems (evidently discarding to fractional derivative and other standard 'procedures' used in fractal Pharmacokinetics), considering that drug release dynamics can be described through continuous but nondifferentiable curves (multifractal curves) is proposed. Then, instead of 'working' with a single variable described by a strict, nondifferentiable function, it is possible to 'operate' only with approximations of these mathematical functions, obtained by averaging them on different scale resolutions. As a consequence, any variable purposed to describe drug release processes will still perform as the limit of a family of mathematical functions, this being nondifferentiable for null scale resolutions and differentiable otherwise.

Preparation of the drug delivery systems
The drug delivery systems were prepared by in situ hydrogelation of chitosan with salicylaldehyde in the presence of DCF, according to a procedure developed in our group (Craciun et al., 2019;Ailincai et al., 2018;Iftime et al., 2020). Using two different molar ratio of the amine/aldehyde functional units, (1.5:1 and 2:1), and keeping constant the amount of drug (1.5 mg), two formulations coded D1.5 and D2 were prepared. Thus D1.5 formulation was prepared by crosslinking of 43.4 mg chitosan with 17.1 mg salicylaldehyde mixed with 1.5 mg DCF, and D2 was prepared by crosslinking 46.71 mg chitosan with 13.79 mg salicylaldehyde mixed with 1.5 mg DCF.

Methods
The corresponding xerogels of the formulations were achieved after they have been lyophilized using a Labconco FreeZone Freeze Dry System equipment, for 24 h at À54 C and 1.512mbar.
The morphology of the formulations was investigated on the corresponding xerogels, using a field emission Scanning Electron Microscope (SEM) EDAX -Quanta 200 at accelerated electron energy of 20 KeV.
The supramolecular ordering of the formulations was investigated on the corresponding xerogels by polarized light microscopy (POM) with a Leica DM 2500 microscope.

The in vitro DCF release protocol
The in vitro drug release has been monitored over 10 days, simulating the in vivo physiological conditions, by using phosphate buffer solution (PBS) of pH 7.4 and keeping a constant temperature at 37 C on the entire investigation period. The formulations of similar amounts (62 mg) were prepared as pills by pressing into a hydraulic press (2 N/m 2 ). The pills were dipped into vials containing 10 mL of PBS. At fixed intervals, 2 mL aliquots were withdrawn and replaced with 2 mL fresh buffer. The supernatant samples were collected and the DCF amount was determined by quantitative absorption spectroscopy, recording the characteristic absorption band at 275 nm, and fitting its absorbance on a predetermined calibration curve (Iftime et al., 2020;Dragan & Cocarta, 2016). The cumulative release of the DCF was estimated from the Lambert-Beer law. The UV-Vis spectra were recorded on an UV-visible spectrophotometer (Perkin Elmer, Lambda 10).

Theoretical consideration
In order to assess the mechanism of the DCF release from the formulations on the two stages, the data were fitted on the following models (Ritger & Peppas, 1987;Masaro & Zhu, 1999;Siepmann & Peppas, 2011): where Q t is the quantity of DCF dissolved in the time t and K 0 is the zero order release constant. First order model: logQ t 5 k · t/2.303, where Q t is the amount of DCF released in the time t and K is the first order release constant. Korsmeyer-Peppas model: M t /M ' 5 k · t n , where M t /M 1 is the fraction of DCF released at the time t, K is the release rate constant and n is the release exponent. Higuchi model: Q t 5 k H · t 1/2 , where Q t is the amount of DCF released in the time t and K H is the Higuchi dissolution constant. Hixson-Crowell model: Wo 1/3 -Wt 1/3 5 k · t, where W 0 is the initial amount of DCF in formulations, W t is the remaining amount of DCF in formulation at time t and K is a constant.

Short reminder on nondifferentiability calibrated on drug release process
Considering that the dynamics of the polymer-drug structural units take place on continuous but nondifferentiable curves (multifractal curves), these dynamics will be described through the scale covariance derivative (for details see Merches & Agop, 2016): In the abovewritten relations, x l is the fractal spatial coordinate, t is the nonfractal time having the role of an affine parameter of the motion curves,V l is the complex velocity, V l D is the differential velocity independent on the scale resolution dt, V l F is the nondifferentiable velocity dependent on the scale resolution, D lp is the constant tensor associated with the differentiablenondifferentiable transition, k l þ k p þ À Á is the constant vector associated with the backward differentiablenondifferentiable drug release k l À k p À À Á is the constant vector associated with the forward differentiablenondifferentiable drug release processes, f a ð Þ is the singularity spectrum of order a of the fractal dimension D F and a is the singularity index. There are many modes, and thus a varied selection of definitions of fractal dimensions: the fractal dimension in the sense of Kolmogorov, the fractal dimension in the sense of Hausdorff -Besikovitch etc. (Mandelbrot, 1982;Jackson, 1993;Cristescu, 2008). Selecting one of these definitions and operating it in the drug release dynamics, the value of the fractal dimension must be constant and arbitrary for the entirety of the compartmental analysis: for example, it is regularly found D F < 2 for drug release correlative processes, D F > 2 for drug release noncorrelative processes, etc. In such a conjecture, through (3) it is possible to identify not only the 'areas' of the drug release dynamics that are characterized by a certain fractal dimension (mono-fractal drug release dynamics), but also the number of 'areas' whose fractal dimensions are situated in an interval of values (multifractal drug release dynamics). More than that, through the singularity spectrum f a ð Þ it is possible to identify classes of universality in the drug release dynamics laws, even when regular or strange attractors have different aspects.
If the drug release dynamics are described through Markovtype stochastic processes (Mandelbrot, 1982;Jackson, 1993;Cristescu, 2008): and for where k is a specific coefficient associated to the fractalnonfractal scale transition and d il is Kronecker's pseudotensor, the scale covariant derivative becomes: In the context for drug release dynamics described by means of motions on Peanotype curves, which implies D F ¼ 2, the scale covariant derivative (5) takes the standard form from the Scale Relativity Theory (Nottale, 2011): where k D is the diffusion coefficient associated to fractal nonfractal scale transition. Therefore, this model, generalizes all the results of Nottale's theory (i.e., Scale Relativity Theory) (Nottale, 2011). Now, accepting the functionality of the scale covariance principle, i.e., applying the operator (1) to the complex velocity fields (2), in the absence of any external constraint, the motion equations of the polymerdrug structural units dynamics (i.e., the geodesics equation on multifractal space) takes the following form: This means that the multifractal acceleration, o tV i , the multifractal convection,V l o lV i and the multifractal dissipation D lk o l o kV i make their balance in every point of any multifractal curve (Irimiciuc et al., 2018) of the polymer drug structural units dynamics. Particularly, for (3) and (4), the motion Equation (7) becomes: Now, separating the polymer-drug structural units dynamics on scale resolutions (the differentiable and nondifferentiable scale resolutions), (7) becomes: while (8) takes the form: For irrotational motions of the polymer-drug structural units, the complex velocity fields (5) take the form: where W is the states function. From here, where ffiffiffi q p is the amplitude and s is the phase, the complex velocity fields (11) become explicit: which enable the definition of the velocity fields: By (12), (14) and (15) and using the mathematical procedures from (Merches & Agop, 2016;Agop & Paun, 2017), the geodesics Equation (13) reduces to the multifractal hydrodynamictype equations: with Q the specific multifractal potential: The Equation (16) corresponds to the specific momentum conservation law of multifractal type, while Equation (17) corresponds to the states density conservation law of multifractal type. The specific multifractal potential (18) implies the specific multifractal force: which is a measure of the multifractality of the motion curves.

Synchronous drug release phenomena
Let it be considered the onedimensional multifractal hydrodynamictype Equations (16-18), in the form: These equations for the initial and boundary conditions: (23) with V 0 the initial velocity and a the parameter of Gaussian distribution of positions, using the mathematical procedures from (Merches & Agop, 2016;Agop & Paun, 2017), admit the solution: the multifractal degree. From here, through (15) the nondifferentiable velocity V F takes the form: Introducing the nondimensional variables: and nondimensional parameters with s 0 the specific time, (24), (25) and (27) become: q n, g, l, / ð Þ¼ p 1 2 aq x, t, r, a ð Þ ¼ 1 þ l 2 g 2 À Á À 1 2 exp À n À g ð Þ 2 Now taking out the quadratic term in g between (30) and (32), it results that for n ¼ const: the ratio V U is homographic dependent of n by the form: From here, the condition (dynamical simultaneity): (i.e., the extension of the first principle of Newton to any scale resolution, or equivalently, 'synchronizations' of drug release dynamics at differentiable scale with drug release dynamics at nondifferentiable scale), implies correlations between phase and amplitude of the shape function, by the form: where q 0 and s 0 are integration constants. Thus, it is stated that various 'mechanisms' involved in the drug release process can be mimed through period doubling, quasiperiodicity, intermittences etc. (for details see Ailincai et al., 2020). Because through the restriction (34) given, for example, by V ¼ ÀU, the multifractal type conservation laws (20) and (21) take the form of the multifractal type 'diffusion' equation: it results that these 'mechanisms' 'manifest'/are 'perceived' as diffusions at various scale resolutions in a multifractal space (fickiantype diffusion, nonfickiantype diffusion etc.) To explain such a situation: the onedimensional drug diffusion of multifractal type from a controlledrelease polymeric system with the form of a plane shut, of thickness d: If drug release of multifractal type occurs under perfect sink condition, the following initial and boundary conditions can be assumed: where q 0 is the initial drug states density of the multifractal type in the 'device' of multifractal type and q 1 is the drug states density at the 'polymerfluid' interface of multifractal type. This solution equation under these conditions can take the following form (for details in the classical case see Crank, 1965). In Figure 1 there are represented the An accurate expression can be obtained for small values of t since the second term of (38) disappears and then it becomes: In such a context, q t q 1 can be assimilated to the fraction of dissolved drug i.e., Mt M 1 q t q 1 , where M t is the amount of drug dissolved in time t and M 1 is the total amount of time dissolved when the pharmaceutical dosage form is exhausted.

Results and discussion
In view of modeling the drug release characteristic, two formulations based on chitosan, salicylaldehyde and DCF (noted D1.5 and D2) were prepared applying the procedure of the in situ hydrogelation described in Section Materials and methods. The formulations have different crosslinking density (NH 2 /CHO ¼ 1.5 or 2, respectively) and the same amount of encapsulated drug. As the formulation impacts the drug delivery by both, size and distribution of the pores of hydrogel on a hand, and size and distribution of the drug crystals on the other hand, the formulation morphology was evaluated by scanning electron microscopy (SEM). As can be seen in Figure 2(a,c) the samples revealed a porous morphology, with interconnected pores, without visible drug crystals, indicating that the DCF was encapsulated into the pore Figure 1. 3 D (left-side) and contour plot (right-side) representations of our multifractal function (Equation (38)) use for drug release mechanism analysis. walls at submicrometric level. This observation was further confirmed by polarized light microscopy (POM) which displayed strong banded birefringence, characteristic for layered ordered phases (Marin et al., 2009) and no crystals of drug into the pores walls (Figure 2(b,d)). This reinforces the idea that the DCF drug is intimately mixed with the chitosanbased matrix, probably by intermolecular forces, forming 'fractal' drug-polymer structural units. Starting from this hypothesis the DCF release can be seen as a progressive release of the fractals of different size.  The in vitro release investigation revealed a progressive release of the DCF in two well delimited stages, as can be seen in the graphical representation depicted in Figure 3(a): a burst release during 8 hours in the first stage, and a slower, continuous release during 10 days in the second stage. The two stages look slightly different for the two samples, reflecting the influence of the crosslinking ratio and the hydrogelation speed. Thus, in the first stage D1.5 formulation released 33% DCF, while D2 only 25%. On the contrary, in the second stage, a more rapid release was noted for the D2 formulation. Both samples reached a total drug release of 70% after 10 days. In the light of the fractal formation hypothesis, it can be appreciated that in the case of D1.5 with a faster hydrogelation rate of the more viscous system, the drug was most probably encapsulated as bigger crystals , forming a larger fraction of fractals richer in DCF, more prone for dissolution in the release medium in the first stage. This can explain their faster dissolution in the first stage, followed by a sustained release over 10 days, of the larger remnant fraction of the fractals scarcer in DCF. In the case of formulation with lower crosslinking ratio, D2, the slower hydrogelation rate allowed a finer distribution of the DCF into the hydrogels matrix forming a larger fraction of scarcer DCF fractals, slowing down the release in the first stage and supporting a sustained release in the second stage. In  Figure 3(b,c) we have compared the experimental date for both the D1.5 and D2 with the predictions made by our multifractal drug release model presented in the previous section. The calibration of the model onto the experimental particularities and the biochemistry of the drugpolymer matrix are performed by fitting the empirical data with the multifractal function. We do notice a good correlation between the theoretical simulated prediction and the empirical data.
The kinetics release mechanism of the DCF drug from D1.5 and D2 samples was assessed by fitting the in vitro release data on the mathematical equations of the Korsmeyer-Peppas, Zero order, First order, Higuchi and Hixson-Crowell, on each of the two stages ( Figure 4, Table 1). As it can be seen in Figure 4(a,b), all in vitro release data proved a good fitting in the first stage (Figure 4(a)) and second stage (Figure 4(b)) on the all five mathematical models. This indicates that the DCF release mechanism is controlled by both dissolution velocity and diffusion through the hydrogel network. This conclusion fits very well with the findings of the theoretical mathematical model, converging to the conclusion of a multifractal composition of the drug release systems which determines a progressive release of the dominant fractal population.
Some correlations can be made between the multifractal and classical approaches; however these need to be understood as mostly qualitatively, considering the early stage of development for our multifractal model. The multifractal model can offer as output parameters fractalization dimension or scale resolution which can potentially be correlated with the transport constants. If we implement the fractal analysis in a sequential manner, as classical models do, we obtain for the first stage a fractalization degree between 2.4 (D2) and 1.6 (D1.5), while for the second stage we get value about one order of magnitude lower 0.7 (D2) and 0.4 (D1.5). This is in line with the differences seen in the constant transports given by Korsmeyer-Peppas or Hixson-Crowell models. Therefore, the transport constants can be corelated with the fractality degree of the drug release process. There is also a clear difference between the fractalization degrees for the D2 and D1.5 which is a direct reflection of the increased number of interactions in the drug release exchange.

Conclusions
Two diclofenac release systems based on a chitosan were prepared by an in situ hydrogelation process with different amounts of salicylaldehyde crosslinker, leading to systems with different crosslinking degree. The in vitro release and fitting on traditional models revealed a prolonged drug delivery mainly controlled by dissolution velocity and diffusion through the matrix. A different drug release rate noticed in the first release stage was attributed to the formation of polymer-drug structural units of different size, controlled by the system viscosity during the hydrogelation step. This indicated the formation of a multifractal system during the hydrogelation which controls the rate of drug release.
A theoretical model in a multi fractal paradigm was developed for understanding the drug release dynamics, considering that these behaviors are described by continuous but nondifferentiable curves. In such a context the irrotational type dynamic of the polymer drug structural units implies the functionality of a multifractal type hydrodynamic formalism. For the unidimensional case of this multifractal type hydrodynamic formalism, it was seen that ratio between the differentiable velocity and the nondifferentiable one for a certain distance depends on a homographic manner on time. The conditions for the simultaneous dynamics imply the synchronization of the drug release mechanisms at the two scale resolutions, expressed through diffusion functions of multifractal type (the diffusion process depends on the scale resolutions).

Disclosure statement
No potential conflict of interest was reported by the author(s).