Comparative statistical analysis of the release kinetics models for nanoprecipitated drug delivery systems based on poly(lactic-co-glycolic acid)

1For more than 50 years, through nanotechnology and biomedicine, formulations and devices have been developed with the aim of improving the delivery and transport of drugs inside the human body. The use of drug delivery systems (DDS) has brought great benefits such as altered drug pharmacokinetics, increased bioavailability, controlled and sustained drug release, target therapy, combinatorial therapy, versatility in presentation modalities, simplified treatment schedules, suit clinical needs and patient preferences, decrease side effects, and increase patient treatment compliance. DDS generally include: microsized forms as polymeric microspheres, microsponges, microneedles, and nanosized devices like liposomes, niosomes, metals nanoparticles, magnetic nanoparticles, nanocrystals, organic nanotubes, exosomes, and polymeric nanoparticles. Especially, polymeric nanomaterials have attracted attention as DDS due to their excellent biocompatibility and biodegradability properties, as well as the spatio-temporal control of drug delivery. Another attractive feature is that they allow the encapsulation of drugs with low water solubility and/or a short in vivo life, thus increasing their bioavailability over time.

2One of the synthetic biodegradable polymers used in DDS is Poly (lactic-co-glycolic acid) or PLGA. This is a biopolymer synthesized from the polymerization of polyglycolic acid (PGA) and polylactic acid (PLA) monomers. PLGA is stable, biocompatible, degradable under physiological conditions, non-immunogenic, and non-toxic, and commercially available with a wide range of different physicochemical properties. The selection of specific characteristics such as its molecular weight, inherent viscosity or the LA:GA ratio allows to control the duration and behavior of the release profile of a PLGA DDS. These particularities make it suitable in applications such as tissue engineering, bioimaging, therapeutic devices, and pharmaceutical formulations as an excipient or as part of DDSs. Therefore, several PLGA-based products are approved by the FDA and the European Medical Agency and are marketed to treat different conditions.

3PLGA-based DDSs include morphologies such as microparticles or microspheres, nanofibers, films, foams, gels, nanoparticles, among others. Being of a size that allows them to navigate freely through the body and remain for a long period of time in the blood circulatory system, nanoparticles, especially polymeric ones, have shown great potential as DDS. For instance, for long-term administration formulations, the time of permanence in blood varies between 24 hours (i.e. Lupron®, AbbVie Inc.), 4 weeks (Lupron Depot®, AbbVie Inc.), 12 weeks (Lupaneta Pack®, AbbVie Inc.), to 24 weeks (i.e. Triptodur®, Arbor Pharmaceutics). To optimally develop a drug delivery system, the release mechanisms and physicochemical processes involved in the transport and release of the active ingredient within an organism must be known, especially when working with nanometric systems. The main role of PLGA in these systems is the regulation of drug release kinetics in order to achieve a sustained and controlled release. This is made possible due to the ability to modify the biodegradability of the polymer by altering its chemistry (lactide/glycolide ratio), molecular weight, stoichiometry, functionalization of the terminal groups, and by changing its size, shape, and porosity. In addition, release disturbances caused by the particle’s ability to absorb water and polymer-drug and drug-drug interactions may occur. Drug transport systems made with PLGA nanoparticles tend to show a biphasic behavior in which initially there is a sudden release dominated by the diffusion phenomenon in which a high concentration of active agent is rapidly released in a short period of time, followed by a period of slow and continuous release, where polymer degradation/erosion is the most influential mechanism. The initial stage is due to the presence of non-encapsulated drug molecules on the surface, drug particles attached to the surface of the nanoparticle, diffusion of the active compound through the polymeric matrix, diffusion of the encapsulated agent through water-filled pores, penetration of the surrounding medium from the surface to the center of the nanoparticle and activation of the hydrolytic degradation of the PLGA. Initial hydrolytic degradation of the polymer and diffusion of the drug through the polymer matrix are also present, although their impact is minor. The hydrophilic or hydrophobic nature of the encapsulating agent and its molecular weight also play an important role in this phase. This initial stage can be shortened or eliminated altogether by adjusting the nanoparticle synthesis technique. The second phase is due to the integration of more complex mechanisms such as the existence of different crystalline and hydrated forms of the drug, the interactions formed between polymer and drug, and the plugging of solution-filled pores, which is the result of the formation of bonds between polymer chains. Sometimes a third and final phase is observed, in which the release increases because the degradation of the polymer within the nanoparticle is catalyzed by the decrease in pH generated by the acidic residues produced by the degradation of the nanoparticle surface. This phenomenon is known as autocatalytic heterogeneous degradation and is frequently observed in large particles.

4Despite the benefits of using micro- and nanoformulations in drug delivery systems, not many products have reached the market. This demonstrates the difficulty of developing PLGA encapsulations and consequently their high manufacturing cost, since at each stage of the design cycle, delivery methods need to be evaluated. The development of advanced DDS depends on a carefully studied and monitored selection of the components. This can be facilitated by the use of mathematical models that allow predicting the rates and behavior of the drug release process, over-time drug concentration in the body, and the biological efficacy of developed DDSs. By choosing the right model it would be possible to reduce the number of experiments required and gain an understanding of the physicochemical dynamic of the phenomena, thus facilitating the development of new pharmaceutical products.

5Mathematical models in general are able to interpret and predict experimental findings and are used as an abstraction of a real system with assumptions and simplifications. These mathematical expressions provide a qualitative and quantitative description of the main phenomena involved in drug release by incorporating more parameters in their equation, related or not to physicochemical properties.

6These models can be: i) mechanistic realistic: real phenomena such as drug diffusion or dissolution, erosion, swelling, precipitation, and degradation of the polymer are taken into account, so they can be complex to apply, or ii) empirical/semi-empirical: not based on real chemical, physical or biological phenomena but are partially realistic from the physical-chemical point of view and are therefore generally easy to apply and explain.

7This study aims to perform a comparative statistical analysis of existing mathematical models of drug release kinetics from PLGA nanoparticle-based DDSs engineered by nanoprecipitation to find which one could be the most suitable general model for describing the release profile.

8Forty data sets describing the over-time drug release of different PLGA nanoformulations as synthesized using nanoprecipitation techniques were extracted from twenty scientific articles. The articles were selected by searching for the following keywords or combinations: “drug release systems”, “poly(lactic-co-glycolic acid)”, “PLGA”, “nanoprecipitation” in academic search engines such as Google Scholar, PubMed y NCBI. To limit the study, the authors selected scientific articles from 2016. Furthermore, to the best of our knowledge, this was the year in which the mathematical model of hyperbolic tangent function was applied for the first time to evaluate the release profile of core-shell lipid nanoparticles. From each data set, the characteristics of the nanoparticle such as encapsulated drug, size, shape, and surface modifications, as well as the characteristics of the polymer used such as LA/GA ratio and molecular weight or inherent viscosity, among others, were listed. Fig 1 illustrates the search process for the scientific articles selected for this study. The data sets obtained can be found in the supplementary information.

9Forty data sets obtained were fitted by linear sum-of-squares regression to five conventional empirical/semi-empirical models:

10Zero-order: The release kinetics of an agent is only a function of time, and the process takes place at a constant rate independent of the concentration of the active agent. Where:

11Qt = percentage of the released drug after time t.

12Q0 = initial percentage of the released drug (Usually 0).

13k0 = rate constant for zero-order kinetics.

14t = time.

15First-order: First-order release kinetics states that the variation in concentration with respect to time depends only on the concentration of the active agent remaining in the device. Where:

16Qt = percentage of the released drug after time t.

17Q0 = initial percentage of the released drug (Usually 0).

18k1 = rate constant for first-order kinetics.

19t = time.

20Korsmeyer-Peppas: It is used to describe drug release from a polymeric system considering non-fickian mechanisms. The model is useful when the release mechanism is unknown or when more than one type of drug release phenomenon is involved. Even if Peppas (1985) points out that this model is applicable only for the first 60% of the release curve; however, we analyzed it in both cases, 60% and 100% drug release profile. Where:

21Qt = percentage of the released drug after time t.

22kKP = constant of nanoparticles incorporating geometric characteristic structures.

23n = release exponent (related to the drug release mechanism).

24t = time.

25Weibull: It is an empirical model that serves as a theoretical basis for almost all release kinetics in heterogeneous matrices. It represents a distribution function describing a phenomenon associated with a finite time. Where:

26Qt = percentage of the released drug after time t.

27α = scale parameter, defines the timescale of the process.

28β = curve shape factor.

29t = time.

30Hyperbolic tangent function: This function arises from the modification of the diffusive release model used by Korsmeyer-Peppas with the intention of being able to extend the mathematical expression to the entire release curve of a homogeneous particle. Where:

31Q∞ = total percentage of the released drug.

32Qt = percentage of the released drug after time t.

33α = constant related to particular size and diffusion constant.

34Fig 2 shows a time-line of the development of the mathematical models used in this study.

35In addition, the 40 data sets were fitted by regression to second to fifth order polynomials, according to the equations shown in S1 File.

36The fitting process of each model was carried out by means of simple linear regression analysis by sum of squares in the free software R. Model acceptance criteria, such as the linear regression indexes (intercept I and slope S), regression coefficient (R2), adjusted or corrected regression coefficient (R2a), sum of squares of residual (SSR), sum of squares of the error (SSE), sum square of total variation (SST), the T-statistic and the F-statistic were calculated using the function “lm” in R.

37The significance test of the t-statistic (t-value) was used to evaluate the capacity of the parameters to linearly relate the dependent variable to the independent variable. In this test the following hypotheses were considered:

38Considering βi as the i-th analyzed parameter. Following the rejection rule, for a significance level of ∝ = 0.05, in a two-tailed test, the critical value according to the probability table is |t| > 1.96, therefore:

39By rejecting H0, it is concluded that βi ≠ 0 and hence a statistically significant relationship between the variables analyzed.

40In case of the intercept, when the null hypothesis was accepted by means of the critical value or the p-value, it was recalculated considering the intercept as 0. If the null hypothesis is accepted for the slope, then the model was rejected.

41On the other hand, the significance test of F-statistic, like t-statistic, indicates whether there is a significant relationship between the independent and dependent variables. In this sense, if only one parameter is analyzed, it will lead to the same conclusion as the t-statistic, but, having more parameters, F-test can prove that there is an overall significant relationship. In this case, same hypotheses were assumed:

42The rejection criteria were:

43By rejecting H0, it is concluded that there is a statistically significant relationship between the analyzed variables, conversely, if p − value > 0.05 the model was rejected.

44Additionally, using the R function “glance” of the package “broom” the value of the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) were calculated for each of the fits.

45The equations for the mathematical model acceptance criteria are summarized in S1 File.

46Whereas the data sets obtained from literature contain few cumulative percentage values over time (3–14 per set), Bootstrap resampling technique was employed. This is a resampling method for statistical inference when the original sample is small. The technique seeks to create multiple data subsets from a real data set without making arbitrary assumptions, with the new data set being equally representative of the target population. It is also commonly used to estimate confidence intervals, bias, and variance of an estimator or to calibrate hypothesis tests, see Fig 3.

47In this study, the Bootstrap resampling technique was used to estimate the values of intercept, slope, and R2 with 1000 Bootstrap samples using the “Boot” function. This requires the packages: “carData”, “car”, and “MUMIn”. To correct the AIC value, we applied the function “boot.stepAIC” of the libraries “bootStepAIC” and “MASS”. In addition, the confidence interval for each of the parameters corrected to 95% was calculated with the function “confint” (S3 File).

48The cumulative drug release data extracted from the scientific articles collected are described in detail with the data in S2 File. Also, Table 1 briefly describes the polymer used for the synthesis, the drug encapsulated, the main characteristics of the nanoparticles obtained, the parameters under which the release assay was performed (release medium and duration) and the percentage of drug release achieved in each case. It can be observed that in 80% of the tests carried out, PLGA with an L/G ratio of 50:50 was used, being the polymer with a molecular weight of 7–17 kD the most commonly used. One hundred percent of the nanoparticles generated are nearly spherical and have a size range between 60 and 550 nm, and only 25% of the formulations performed surface modifications. On the other hand, the drug release conditions were variable. 32 of the release assays were performed in Phosphate buffer saline (PBS), 3 in hydrochloric acid, 2 in Phosphate buffer (PB), 2 in Fetal bovine serum (FBS) and 1 in Simulated body fluids (STF). The time to study drug release ranged from 3 hours to 21 days (504 hours). Finally, the number of observations per set was also variable. 20% of the data presented 6 observations, 17.5% 8 observations, 15% 14 observations, 12.5% 7 observations, 10% 5 observations, another 10% 9 observations, 5% 3 observations and only one data set presented 4, 10, 11 and 12 observations respectively.

49S3 File shows, by means of an example using set number 33 and its fit to the Korsmeyer-Peppas model, the program used to fit the mathematical models in the statistical program R.

50S4 File shows the results of the fit and the acceptance criteria for each studied mathematical model and polynomial. The first criteria for rejecting a mathematical model was to evaluate the slope and F-statistic, after recalculating with intercept = 0 if necessary. Table 2 shows the results obtained. In the case of zero-order fit, sets 7–10, 18 and 19 were rejected, for first-order sets 5–11, 18, 19, 28, 29, 34 and 40, for Korsmeyer-Peppas 100% sets 7–9, 18 and 19, for Korsmeyer-Peppas 60% sets 7–10, 18, 26, 27, 28 and 36, for Weibull sets 7–9 and 18, and for Hyperbolic Tangent only the fit for sets 7 to 9 was rejected. It is important to note that no model was adequately fit for sets 7, 8 and 9. In some cases the number of observations limited the application of the statistical analysis in R. Specifically for the Korsmeyer-Peppas 60% model, it was not possible to perform the analysis on sets 7–10, 18, 26, 26, 27 and 36 because with the two first observations the set exceeded 60% of the release curve. We analyzed R2- R2a values, and AIC-BIC of previously un-rejected models. The best fit for each dataset is detailed in Table 3. According to these criteria, none of the datasets fit the zero-order and first-order models. Neither model fitted Weibull according to AIC and BIC. When studying in parallel R2- R2a and AIC-BIC only 48% of the sets accept the same model under these two criteria, which indicates that they are not equivalent. The mathematical model with the largest number of data sets with adequate fit is the hyperbolic tangent. Following R2- R2a and AIC-BIC criteria, 48% and 65% of the sets fit this model best, respectively. In addition, sets 7, 8 and 9 sets (8% of the data) did not fit any model.

51Regarding the fit to the polynomial equations, the value of the intercept was analyzed (a letter in all polynomial functions) following the criteria mentioned in the methodology section. If null hypothesis was accepted, then the statistics were recalculated with a value of intercept = 0. If the value of the last term, represented by the letters c, d, e, f for each polynomial degree respectively, was equal to zero the model was rejected since the equation is transformed into one of a lower degree. The results are described in Table 4. Through this analysis it was found that in this particular case, the results obtained by AIC-BIC y R2—R2a are the same for all sets. Nine of the forty datasets did not fit any model. In accordance with the t-value and p-value of the last term, 11 data sets were rejected for the second-degree equation, 23 for the third-order equation, 27 for the fourth-order equation, and 25 for the fifth-order equation. AIC-BIC and R2—R2a showed that the polynomial equation that provided the best fit to the largest amount of data was the second-degree equation (38%), followed by the third-degree equation (20%), the fifth-degree equation (13%) and finally the fourth-degree equation (8%). It is worth noting that when the average R2—R2a of the data sets unrejected previously were evaluated, it was found that these values tend to increase as the degree of the polynomial increases, as can be seen in Fig 4.

52Moreover, when Bootstrap resampling was performed (S5 File) it could be found that the best fit for data sets 7, 8 and 9 is Korsmeyer-Peppas 100% or zero-order by R2 and Korsmeyer-Peppas 100% by AIC, see Table 5. More specifically, by AIC, 30% of the sets present a better fit with Korsmeyer-Peppas 100% and 70% with Hyperbolic Tangent for polynomial equations 13%, 23%, 23%, 23%, and 33% were fit to second to fifth degree polynomials respectively. In addition, a single data set fitted to all polynomials and 3 to none. It is important to note that by analyzing the coefficient of determination of the resampled data, more than one model is accepted. As the best fit, by using the AIC with the most negative value, the model or polynomial equation that best fits the data can be determined.

53Finally, Figs 5–7 show scatter diagrams of R2 in relation to AIC and BIC for the models analyzed, both for the values obtained by linear regression and the Bootstrap corrected values. Additionally, the region of what is considered a good fit is indicated: between 0.8 and 1 for R2 and between 0 and -inf for AIC and BIC. In both figures it can be clearly observed how in the Hyperbolic Tangent Function model the data sets are mostly grouped in the region considered as having the best fit.

54In this study, five general controlled drug release models and four polynomial equations, from the second to the fifth degree, were statistically analyzed. The main objective was to find a mathematical expression that can describe the release mechanism of an active compound observed in PLGA nanoparticles synthesized by nanoprecipitation method.

55While it is true that mathematical models are tools that help to understand the dynamics of a drug release, and to save time and money, it is important to note that these also have shortcomings in terms of predictive power and in explaining and addressing all the physical, chemical and biological phenomena involved since they mainly describe the dominant event. Among them is the initial burst release, as the fraction released can sometimes approach 80% of the total encapsulated agent in a time that may be less than 10% of the duration of the entire release process. In addition, very few models described in the literature consider the effect of drug properties, especially acids, on the degradation rate of polymer particles. It should not be forgotten that mathematical models are approximations of reality where assumptions are made, so they should always be accompanied by experimental data that can support these assertions. Empirical and semi-empirical mathematical models tend to have limited predictive power, so it is recommended that these theories be applied, and be of great use, when you want to compare release profiles using a specific parameter.

56On the other hand, the coefficient of determination “R2” is one of the most widely used criteria to determine whether or not a model is a better fit to a release curve. However, it has been widely demonstrated that R2 can increase, i.e., tends more to one, if more variables are added to the equation, even if they are completely irrelevant to the model. This problem is more accentuated when the number of observations is insufficient. Similarly, R2 is not considered an appropriate goodness-of-fit criterion when analyzing multi-parameter models, as our study. Therefore, taking into account the adjusted coefficient of determination R2a is the most appropriate, although in a limited way, since its value only increases when a relevant variable is added to the model, otherwise it decreases. Since there is no perfect correlation between the variables, special care should be taken in the interpretation of the results when R2 or R2a take the value of 1. In addition, it is recommended to accompany these criteria with statistical procedures such as t-test, F-test, and prediction intervals to know the quality and influence of individual observations in the estimation.

57Since the Bayesian and Akaike criteria are based on the addition of statistical errors corrected by a penalty function proportional to the number of parameters, they can be a good alternative to determine which model offers the best fit. In addition, these criteria help to determine whether the inclusion of a term, or increased complexity in the equation, results in better fit. It is important to emphasize that models with more negative values of AIC and BIC are likely to minimize information loss, which translates into a more appropriate choice of the model with the best fit.

58The zero-order model is usually employed to assume a single-phase ideal condition in which the release kinetics is continuous and can be explained as a function dependent on time and independent of the dissolved concentration in the release medium. In addition, zero-order behavior can be observed in materials capable of swelling, yet due to complicated bulk erosion and surface erosion and diffusion processes, most studies on the sustained release of hydrophilic drugs have failed to achieve a release profile that conforms to zero-order kinetics for a prolonged time. This conclusion was verified in this study, in which the data sets had a biphasic behavior far from zero-order according to the low values of R2 and R2a and the high averages of the calculated AIC and BIC. Likewise, the results obtained showed rejection of the first-order model for all data sets, demonstrating that the data sets do not follow a model that is dependent on the initial drug concentration.

59The Korsmeyer-Peppas equation also called "power law" is a semi-empirical model based on the diffusion phenomenon that is used to describe in a general way the main transport phenomena involved in the release, either by diffusion or swelling. Although a good fit has been achieved with this model, several authors have shown that Korsmeyer-Peppas is applicable only to the first 60% of the release profile and is not an appropriate equation to describe the complete release profile. By examining only the first 60% of the data sets, some problems could be detected, such as the small number of observations in some sets, which made it difficult and even impossible to fit the data to the models. It was found that 4 data sets presented a more negative AIC value when 60% of the curve was analyzed, and 6 when 100% was analyzed. This demonstrates that although the model is more complex than the previous ones, this does not result in a better description of the release process.

60The empirical Weibull model is an expression that tries to find a linear relationship between the logarithm of drug release and the logarithm of time. Despite not having a kinetic nature, nor parameters physically related to the release phenomena, the ability of the logarithmic form to dampen the abrupt changes in drug concentration rate over time as they pass from one phase to another, has allowed this model to be successfully applied in several release studies. In this study, the Weibull equation was able to fit 12 data sets adequately according to the R2 and R2 adjusted, however, when studying the value of AIC and BIC, the model failed to have the best fit in any of the sets, probably because as an empirical model, its predictive power is poor. Weibull has two geometric parameters. The first one is the scale factor "α" which is part of the intercept with the vertical axis, expressed as "ln(α)". The second one is the slope of this equation, usually called "β", is a shape factor of the release curve. Thus, if β = 1, the curve maintains an exponential shape, if β<1 it is considered parabolic and in the case of having a β>1, the shape of the curve is sigmoidal. When analyzing this parameter in the collected data sets (S6 File), a tendency towards a parabolic shape could be observed by obtaining an average of β = 0.71. Sets 8, 9, 18,19, 33, 38 and 39 were the only ones with a sigmoidal curve (β>1). In order to relate the Weibull model to the predominant transport mechanism in drug release, relationships between the Weibull parameter “β” and the Korsmeyer-Peppas parameter “n” have been proposed. Following this theory, Fickian diffusion would be demonstrated when β≤0.75, anomalous transport (diffusion and swelling) when β is between 0.75 and 1, case II type transport when β = 1 and case II type transport which is the combination of phenomena such as diffusion, erosion and macromolecular ratio of polymer chains for the case of cylindrical tablets with β>1. Doing this analysis with "n" and "β" obtained in this study (S6 File) we obtain a value of “β” to check Fickian diffusion of 0.75, while the transport type case II is explained by β = 1.45. The differences in the ranges could be explained by the difference in the geometry of the DDS. This analysis together with the Korsmeyer-Peppas "n" parameter tests the influence of Fickian diffusion on drug release from PLGA nanoparticles.

61The hyperbolic tangent function model was the one with the best fit to the experimental data analyzed. This mathematical expression was proposed by Eltayeb et al. in 2016 and is derived from the Korsmeyer-Peppas model, considering diffusion as the main transport phenomenon involved in drug release from nanoparticles. This equation uses the hyperbolic tangent of time as a way to mitigate the abrupt effect of the released drug concentration when its phase changes, similar to the effect of the logarithm in the Weibull model. Unlike Korsmeyer-Peppas, this model can be fitted to the complete release profile and not just to its first 60%. Besides, it is applicable to homogeneous matrices, where the ratio of drug and polymer is equal. From the data sets analyzed, 18 sets were fitted to the hyperbolic tangent model by R2 and R2a and 27 sets considering the AIC and BIC, indicating the model’s ability to describe and predict the drug release profile from PLGA nanoparticles synthesized by nanoprecipitation. Being a relatively new model, the physicochemical significance of its parameters is still not fully understood, however, its authors consider that they are related to particle size and diffusion constant. Furthermore, to our knowledge, this model has not been used by other research groups, however, through the analysis performed in this project it is demonstrated that it can be a useful and statistically proven tool to describe release profiles.

62As for the four polynomial equations, the importance of the statistical tests to reject the models was evident. Although they showed a good fit measured by R2 and R2a, the increase of these does not indicate a good fit of the model since they tend to increase with the number of variables analyzed in the equation, as shown in Fig 2. Also, when working with polynomial models, we always try to minimize the number of parameters used in the model since, if the number of parameters is equal to n-1, where n is the number of points to be analyzed, we would obtain an error of fit equal to zero, which is described as an overfitting. Equally, the values of AIC and BIC were more positive compared to the mathematical models. Since AIC and BIC are a measure of goodness-of-fit conditional on the number of parameters used in the model, this result is congruent, and indicates the inability of the expressions to predict the future behavior of the system. Pourtalebi Jahromi and collaborators (2020) used the second degree polynomial equation in their analysis of the release kinetics from nanoparticles, obtaining results similar to those of this study. Nevertheless, no literature was found in which the third-, fourth- or fifth-degree equations were used as a mathematical model in drug release.

63On the other hand, the number of observations per set limited in some cases the application of the mathematical models and their statistical tests, especially in sets 7, 8 and 9. In these sets the authors reported the release of 100% of the encapsulated drug in 3 observations. According to Jenkins and Quintana-Ascencia (2020), having a minimum of 8 observations would allow significant statistical conclusions to be reached; however, to increase reproducibility, it is recommended that a minimum of 25 observations be used. From the data sets collected, 50% of them have 8 or more observations, with sets number 20–25 having the most number of observations (14 each set).

64When Bootstrap resampling was performed, it was possible to determine the model that best fit the data sets 7, 8 and 9. It proves that this technique can be a powerful tool when there is a limited number of observations. The results with this method were similar to the raw data, obtaining that 70% of the data sets fit the Hyperbolic Tangent model.

65This study demonstrated that, although they are tools that help to understand drug release dynamics, there is no general empirical/semiempirical mathematical expression that can describe, in all cases, the release profile of drugs encapsulated in PLGA nanoparticles synthesized by nanoprecipitation methods. In addition, it was found that although R2 and R2a are the most commonly used criteria to determine whether or not a model fits the data obtained in an experiment, they are not the most appropriate. It was also shown that the Akaike and Bayesian criteria can better reflect the fit of a model since the results are not influenced by the inclusion of new terms or an increase in the complexity of the equation. It was also revealed that the number of observations per set is a limiting factor for the application of different models and the subsequent statistical analysis. Therefore, the Bootstrap resampling technique becomes a very useful technique to solve this drawback. In the specific case of this study, 50% of the sets studied do not meet the requirement of at least 8 observations for the construction of the release curve, therefore these results are not very reproducible and have low statistical significance.

66Furthermore, the analysis employed in this project provided significant statistical evidence to consider the Hyperbolic Tangent Function model as the most adequate and general model to describe the drug release kinetics. This model, unlike the mathematical expression of Korsmeyer-Peppas, could be adjusted to the complete release curve. Since it has kinetic parameters, it acquires greater predictive power than the Weibull model. However, a more exhaustive study of this model is required in order to understand the chemistry, physics, and biology behind it.