Mathematical Modeling for Pharmacokinetic Predictions from Controlled Drug Release Nano Systems : A Comparative Parametric Study

In the present work, several mathematical models well-known in the literature for simulating drug release kinetics are compared using available experimental data sets obtained in real systems with different drugs and nano-sized carriers. Herein, the ÷2 minimization method, is employedconcluding that the Korsmeyer-Peppas modelprovides the best-ûtin all cases. Hence, (i) better understanding of the exact mass transport mechanism(s) involved in drug(s) release, and (ii) quantitative prediction of the drug release kinetics, can be computed.

Nowadays, pharmaceutical industries and registration authorities focus on drug dissolution and/or pharmacokinetic release studies.Mathematical modeling aids at predicting drug release rates, and thus helping researchers to develop highly eûective drug formulations and more accurate dosing regimens saving time and money 1 .Kinetic models describe the amount of drug dissolved "C" fromsoliddosage form as a function of time t, or f = C(t).Since in practice the underlinemechanismis usually unknown, some semi-empirical equations, based on elementaryfunctions (polynomials, exponentials etc), are introduced.Up to now, a signiûcant number of mathematicalmodels have been introduced in the literature [1][2][3] , and in principle, onecan opt to use any of these.So, the question naturally arising herein is: which mathematical model is the best-fit to use for a given nano-system?
In the present work, weattempt toreaddress precisely this question by systematically comparing various existingmathematical models.Already in 2 , it is mentioned that statistical methods can be usedto select a model, and one common method is based on minimization of the co-eûcient of determination R 2 , or if models with diûerent number of parameters are to be compared the adjusted coeûcient of determination R 2 adjusted = 1"(1"R 2 )(N-1)/(N "m) is preferred, where N is the number of experimental points and m is the number of free parameters ofa given mathematical model.
Herein, however, and to the best of our knowledge, it is the ûrstattempt in which themathematical model comparison is done explicitly using concrete experimental data thatcorrespond to diûerent drugs and diûerent nanoparticles; a more realistic approach perhaps.Furthermore, we employed the ÷ 2 minimization method instead of the R 2 coeûcient of determination, resulting indiûerent conclusions as we shall discuss in more detail later on.Thereby, the work is organized as follows:We first presentthe models to be compared as well as the data sets we have used for the analysis.Then,we perform the comparison and present findings and conclusions.A narrative format was deemed suitable for added clarity.

Mathematical Models and Data Sets
We compared the following mathematical 6 renowned models [1][2][3] : On the other hand, the obtained data setsare summarized in the tables below: Tables 1 and 2 relate to a multidrugloaded nanoplatform composed of Layer-by-layer (LbL)-engineered nanoparticles (NPs) achieved via the sequential deposition of poly-l-lysine (PLL) and poly(ethylene glycol)-block-poly(l-aspartic acid) (PEG-b-PLD) on liposomal nanoparticles (LbL-LNPs).The multilayered NPs (<"240nm in size, illustrated in Figure 1) were designed for the systemic administrationof doxorubicin (DOX -release kinetic profiling is displayed in Figure 2) and mitoxantrone (MTX).Data sets in Tables 3 and 4 relate to poly(D,L-lactide-co-glycolide) (PLGA-based nanoparticles) designed for the longterm sustained and controlled (linear) delivery of simvastatin (SMV).Finally, [poly(å-caprolactone)based nanocapsules were prepared for the data set summarized in Table 5.

Model Comparison
We now proceed to perform the model comparison using the c 2 minimization method.For a given data set with N number of time points with values Q i and errors ó i , i taking values from one to N, and for a given function f(t;a 1 ,a 2 ,...,a m ) that models the amount of drug as a function of time and is characterized by m free parameters (where N > m), we compute c 2 using the standard formula: ... (7)  where we sum over all experimental time points from i=1 to i=N, and thus c 2 is a function of the free parameters that characterize the mathematical model.Minimizing This last step is necessary in order to compare models with diûerent number of free parameters.
In our analysis the models are characterized either by one or by twofree parameters, and so m = 1 or m = 2, while the data sets have either 8, 10 or 12 pointsand so N = 8, N = 10 or N = 12.
For a given data set the model that best ûts the data is the one with the lowest ÷ 2 min /d.o.f.We start with the ûrst data set seen in Table 1 and we minimize ÷ 2 for all models one by one using the computer software Mathematica 11 .By comparing ÷ 2 min /d.o.f we see that the power law model has the best ût.The values of the parametersare summarized in Table 6, while as was illustrated in Figure 2, we can see that indeed the power law model fits the data way better than the Higuchi model.We then follow exactly the same procedure for the rest of the data sets seen in Tables2, 3, 4 and 5. Our results show that the power law model has the best ûtin all cases, and therefore our conclusion is robust.
O u r r e s u l t s a r e i n t e r e s t i n g f o r threereasons:Foremost, we have shown that although the most-widely used model in the literatureis the one introduced by Higuchi 4 , at least the class of systems considered here are bestdescribed by the power law model.In addition, we have shown that it is possible that amodel with more parameters has a better ût to the data contrary to what is stated in the literature whenthe coeûcient of determination R 2 is used 2 .This is due to the fact thatalthough the number of degrees of freedom decreases when the number of free parametersincreases, in some cases the c 2 at the minimum is reduced so much that overall the c 2 /d.o.fis lower.Finally, knowing the model that describes the systems studied here in, it would  1. Shown are the data points, the Higuchi model (red color) and the power law model (black color) which fits the data better than the Higuchi model beinteresting to try to understand the underlying mechanism starting from basic principles,and relate the parameters of the model with properties of the system.In that case, sincethe parametersof the model have been already determined upon comparison with the data,one can compute the properties of the system, and thus the properties of the system couldbe measured experimentally using our method.Furthermore, it is interesting to note at this point that the power law time dependencecan be mathematically derived as the exact analytical solution of the diûusion equationin one dimension in the semi-inûnite domain x > 0: C(t,x) t = D C(t,x) xx ...( 8) w h e r e t h e s u b i n d e x t d e n o t e s differentiation with respect to time, while the subindex xx denotes double differentiation with respect to space, with the initial condition C(t = 0,x) = 0 and boundary condition C(t,x = 0) = kt n/2 .In the above initial/boundary problem D is the diûusioncoeûcient assumed to be a constant, C(t,x) is the drug concentration asa function of time andposition and k,n are constants.It is known from mathematicalphysics that this boundary/initialvalue problem is well posed and it has a unique solution 11 .Using the method of Laplacetransform (see e.g. 11) one ûnds that the unique solution that satisûes thediûusion equation and all conditions is the following 12 : C(t,x) = kΓ(1 + n/2)(4t) n/2 i n erfc(x/2√Dt) ...( 9) where Γ(z) is the Euler's Gamma function, and we make use of the error function erf(x)and the complementary error function erfc(x) deûned as follows: ... (10)   ...(11) For more details on the special functions of mathematical physics see e.g. 13 .Finally, given the drug concentration, we can now compute the amount of the drug as a function of timeby performing the integral over all space from zero to infinity: ... (12)  The integral can be computed exactly and ûnally we obtain:

CONCLUSIONS
In this work, we conducted comparisons between several mathematical models widelymentioned in the literature regarding predicting overall release behavior.We have used 5 diûerent data setsobtained experimentally in realistic systems with diûerent drugs and nanoparticles.
Eachmodel is characterized by one or two free parameters to be determined uponcomparison with the data.We have used the c 2 minimization method to determine the values of the parameters of each model, and we have obtained the minimum value of c 2 per degree offreedom for each model.Our results show that among all mathematical models studied herein, the power law model has the best-ût in all 4 cases.We conclude that at least the class ofsystems considered here are best described by the power law model, characterized by two free parameters, although the Higuchi model is the most widely-used in the literature, and also despite other claims that adopting the coeûcient of determination R 2 , models with more parameters have a worse ût to the data.Finally, our derived method could in principle be used to measure variable properties of the nano-systems, experimentally.
c 2 we determine the values of the parameters for which the model best ûts the data, and ûnallywe compute c min 2 /d.o.f,where d.o.f stands for the number of degrees of freedom given by N "m.