Mathematical Modeling of Hydrogels for Drug Delivery Systems PDF

International Journal of Pharmaceutics 560 (2019) 175–190 Contents lists available at ScienceDirect International Journal of Pharmaceutics...

International Journal of Pharmaceutics 560 (2019) 175–190 Contents lists available at ScienceDirect International Journal of Pharmaceutics journal homepage: www.elsevier.com/locate/ijpharm An overview on the mathematical modeling of hydrogels’ behavior for drug T delivery systems Diego Caccavo Department of Industrial Engineering, University of Salerno, 84084 Fisciano, SA, Italy ARTICLE INFO ABSTRACT Keywords: Hydrogels-based systems (HBSs) for drug delivery are nowadays extensively used and the interest in modeling Hydrogels their behavior is dramatically increasing. In this review a critical overview on the modeling approaches is given, Drug delivery quantitatively and qualitatively analyzing the publications on the subject, the trend of the publications per year Modeling and the type of modeling approaches. It was found that, despite the drug release fitting models (i.e. Higuchi’s equation) are the most abundant, their use for HBSs is decreasing in the last years and luckily, considering the limiting assumption on which they were built, they will be confined to simple mathematical fitting equations. Within the mechanistic models the “multi-component” with the swelling approximation (mass transport only) and with the mechanics (fully coupled) are experiencing the highest growth rate, with much more interest toward the last one that, in the next years could be able to provide a first principles model. Statistical models, especially based on the response surface methodology, are rapidly spreading in the scientific community mainly thanks to their ability to be predictive, regardless of the phenomenology, in the analyzed design space with very low efforts. Neural Networks models for HBSs, in countertrend with their use in the pharmaceutical industry, have never take off preferring less data demanding statistical models. 1. Introduction et al., 2000) but, surprisingly, exhibiting a rheological behavior more similar to that of solids (elastic behavior) than to that of liquids (viscous Hydrogels-based systems (HBSs) for controlled drug delivery are behavior) (Lapasin and Pricl, 1995). For all these reasons, gels resemble widely used devices (Caló and Khutoryanskiy, 2015) thanks to the high living tissues much better than any other synthetic bio-material (Lin patients’ compliance (patients prefer oral and/or topical administration), and Metters, 2006). A distinction has to be done between hydrosols and and thanks to the possibility to modulate the active ingredient release hydrogels: with the first are indicated solutions in which polymers are rate, maximizing the therapeutic effects. Modeling the behavior of HBSs dissolved in water, instead, with the term hydrogels cross-linked hy- is currently a topic of great interest since, as in all fields, mathematical drosols are identified (Siepmann et al., 2011). These last, due to the models constitute an indispensable step to reach a deep understanding of presence of chemical or physical cross-links (junctions, tie-points, en- the process, which is necessary for the correct design of these systems, as tanglements) in the polymeric network, are unable to dissolve in water well as for cost reduction (i.e. faster and cheaper experimental cam- but they can absorb it increasing their volume. The behavior of these paign). The wide spreading of hydrogels and HBSs, across several sectors systems is the result of the interaction of all its parts (polymeric net- (from the biomedical to the agro-food to the construction industry) has work, water, ions) with the external environment. Indeed, the water led several researchers, with different backgrounds, to deal with hy- adsorption due to diluting force (the entropy of the system increases, drogel-based systems, producing a vast and fragmented literature on the like in the normal solubilization process of linear polymers) causes the possible mathematical modeling approaches. With these premises, it network swelling. Therefore, the chains between network junctions are appears useful to have a brief overview on the phenomenology of HBSs required to assume elongated configurations generating an elastic force. before starting with the mathematical treatise. As swelling proceeds, this force increases and the diluting force de- creases. Ultimately, a state of equilibrium swelling is reached in which 1.1. Phenomenology these two forces are in balance. When these polymers have an ionic network with ionizable groups, the swelling ability may be greatly in- Hydrogels are three-dimensional, hydrophilic, polymeric networks creased as a result of the localization of charges that, setting up an capable of imbibing large amounts of water or biological fluids (Peppas electrostatic repulsion, tend to expand the network. However, the fixed E-mail address: [email protected]. https://doi.org/10.1016/j.ijpharm.2019.01.076 Received 4 December 2018; Received in revised form 28 January 2019; Accepted 31 January 2019 Available online 11 February 2019 0378-5173/ © 2019 Elsevier B.V. All rights reserved. D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 charges are not the only ions present in the gel, at least a stoichiometric It should be clear by now that hydrogels (or HBSs), when solicited amount of mobile counterions has to be considered. These, by screening by chemical or mechanical stimuli, and therefore forced to a condition effect, reduce the ideal swelling capacity due to the repulsion (Flory, of non-equilibrium, can undergo substantial structure modification 1953). This peculiar behavior, which has led to define this soft matter and/or chemical composition variation to re-establish the equilibrium. as “smart materials”, makes hydrogels and hydrogel-based systems very The mutual interaction between the mechanical structure and the sol- attractive by several frontier fields, such as biomedical applications and vent, apart guiding the driving forces for the mass transport, influence specifically drug delivery applications. In such uses hydrogels are also the transport coefficients. I.e. the polymeric chains exert a sieving loaded with an active ingredient (drug) and eventually with excipients effect on the diffusing molecules, which results in non-constant diffu- at the dry state (mixing and tableting the powders (Barba et al., 2009)) sion coefficients. Generally, the models useful to predict the diffusivity or in the wet state (in solution before the crosslinking (Chen et al., of a molecule “i ” in a swollen network “s ,” Di.s , with respect to the 2004) or by letting the drug diffuse inside the hydrogel (Lin and diffusivity of the same molecule in the solvent “1,” have the following Metters, 2006)). It is clear that in these cases the phenomenology is general form (Lin and Metters, 2006) (the ratio Di.s /Di.1 is sometimes further complicated by the presence of multiple species that can in- called “the retardation effect”): teract and modify the HBSs behavior. In such a complex system, the Di.s mechanics of the structure influences the water absorption/desorption = f( , 2 , rs) Di.1 (1.2) and vice versa, resulting in different mass transport regimes (for the solvent as well as for the other species). These regimes are reported in where, and 2 are, respectively, the network mesh size (is a measure literature as “Fickian” and “non-Fickian” (or “Anomalous”) diffusions. of the space available between the macromolecular chains) and the The first transport is characteristic of a concentration gradient driven polymer volume fraction, while the parameter rs is the size of the dif- processes, whereas the latter deviate from a standard diffusion problem fusing molecule. The mechanistic theories which are used to build the due to the influence of the polymer network relaxation. Vrentas et al. left-hand side of Eq. (1.2) are known as hydrodynamic theories, ob- (Vrentas et al., 1975; Vrentas and Duda, 1977) have suggested a simple struction theories, and theories based on free volume. All these theories way to establish the regions in which Fickian and non-Fickian trans- have been thoroughly reviewed in the following references (Amsden, ports take place, on the basis of the diffusional Deborah number, NDe. D , 1998a,b; Masaro and Zhu, 1999; Grassi et al., 2006). defined according to Eq. (1.1), in which is the characteristic stress- Summarizing, when a dry polymer forming hydrogel, loaded with relaxation time of the polymer–solvent system and D is the character- an active ingredient (AI), is immersed in a physiological fluid (dis- istic time for the diffusion of the solvent in the polymer: solution medium), the solvent penetrates inside the polymeric matrix. If the polymer shows a glass transition temperature higher than the room/ sG(s)ds/ G(s)ds NDe.D = = 0 0 physiological temperature, when the solvent concentration exceeds a D L2Ch/D1.s (1.1) threshold value, polymeric chains unfold so that the glass-rubbery The stress-relaxation time, , can be evaluated by integrals of the transition occurs (at temperatures lower than the pure polymer glass shear relaxation modulus, G(t) , over the entire relaxation time spec- transition temperature) and a gel-like layer, surrounding the matrix dry trum; the diffusion time, D , is given by the ratio between the second core appears (Grassi et al., 2006). The moving front at which this power of a characteristic diffusional path length (for the solvent), LCh , process takes place is called “swelling front”, which separates the and the diffusion coefficient of the solvent in the swollen network, D1.s swollen from non-swollen matrix (Colombo et al., 1999). In the swollen (Davidson and Peppas, 1986). If the change in solvent concentration region the polymeric chains assume an elongated configuration that during the swelling process is limited, average values for each char- allows the contained AI molecules to easily diffuse toward the outer acteristic time have to be used, and then the full process can be char- dissolution medium, once that they are dissolved. Indeed depending on acterized by a single value of the diffusional Deborah number. If the the drug solubility, in the swollen layer there could be zones in which change in solvent concentration is large, the diffusional Deborah the drug coexist in the dissolved and dispersed form (Grassi et al., number have to be calculated for both the initial and final stages, and 2010). The front that separates the swollen matrix, containing only their order of magnitude will be used to characterize the behavior of the dissolved drug, from the swollen part that contains both dissolved and system. The value of the diffusional Deborah number, NDe.D, dis- dispersed drug, is called “diffusion front”. Additionally, on the zone at criminates the nature of the diffusive phenomena: which the swollen matrix is in contact with the outer medium a third front can be defined: the “erosion front”. On this boundary the polymer Large values of the Deborah number (N De.D 1) indicates that the network became extremely hydrated and process like chains disen- tanglement can take place, “eroding” the matrix (Siepmann and characteristic relaxation time, , is long with respect to character- istic diffusion time, D : the polymer structure does not change during Siepmann, 2008). the water diffusion process, i.e., the polymer remains in its glassy state. The diffusion phenomenon is usually described by the con- 1.2. Modeling approaches ventional Fick’s law, using coefficients of diffusion constant and independent from water/polymer concentrations. A first question to answer when approaching hydrogels and HBSs is: Small values of the Deborah number (NDe.D 1) indicates that the “are they multiphasic or monophasic systems”? The answer cannot be taken for granted. Despite in most experimental works the response is relaxation phenomenon is much faster than the diffusion phenom- enon. Practically, it is a diffusion through a viscous mixture (the simply avoided, it become fundamental when the aim is to develop a swollen, rubbery hydrogel), a process which can be described again mechanistic mathematical model of the system. The most natural ap- by the conventional Fick’s law, using coefficients of diffusion which proach is to consider hydrogels as single-phase matter, in which several are a strong function of water/polymer concentrations. components can coexist, like it would be indisputably done for poly- Intermediate values of the Deborah number (NDe.D 1) indicates meric solutions (hydrosols). Another vision is to consider hydrogels as made of different phases, i.e. the water phase (free water) is separated that the two characteristic times are of the same order of magnitude, i.e., the relaxation and the diffusion phenomena take place on the from the polymeric phase (polymer plus bound water), and these can same time scale. This is the transition zone in which the polymer exchange momentum. This approach could sound more representative experiences its glass–rubber phase change, the mixture has a vis- of hydrogels with macro-porosity, whereas the first could appear more coelastic nature, and the diffusion is an anomalous transport (and its suitable of hydrogels with nano/micro-porosity. Mathematically both limit, the Case-II transport), which is non-Fickian. the approaches can be applied, independently on the “real” physical state (porosity) of the system, just choosing a proper representative 176 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 elementary volume (REV) (Caccavo et al., 2017b). As it will be shown 2014), which expensively relies on statistical techniques (i.e. design of in the Section 3, when used for drug delivery purposes, the monophasic experiments). Neural networks models, useful to relate output to input approach is predominant. Another important question is related to the through complex transfer functions (obtained from the “training” of need of modeling/analyze the full behavior, mass transport plus me- layers of neurons), have not been used much in this field and in the last chanics, or just one aspect, mass transport only. The difficulties related years a decreasing trend is recorded and predicted. Reviews on the to the solution/analysis of the full HBS behavior have led many re- “modeling of the drug release from hydrogels” have almost a constant searchers to describe hydrogel-based systems with a “mass transport trend (∼1 each 3 years). only” or, even easier, with empirical equations. From this analysis it can be seen that the topic is alive in the sci- Querying one of the biggest scientific work database (Scopus.com) entific community, predicting a further uptrend in the number of for the articles related to “modeling of the drug release from hydrogels” publications in the next years. Moreover it can be seen that the interest (the query used is in the footnotes1), 337 works (in October 2018) is moving away from the simple fitting equation models, in favor of the appear. Despite the use of “negation operators” in the query string, it more reliable and comprehensive mechanistic and statistical ap- was required a manual screening to eliminate the manuscripts not proaches. strictly dealing with the subject and the ones focusing on the behavior of pure hydrogel (polymer – water only) systems. This operation re- 2. Drug release fitting models duced the database to 153 works, in the year range 1980–2018. Following the previous considerations, these works can be cate- As it has been previously shown, the most relevant part (in number) gorized in: of publications inherent to the “modeling of the drug release from hy- drogels” is based on this approach. In particular, considering the works drug release fitting models; in Table 1, the most used equations can be individuated: mechanistic models; statistical models; the zero order equation; neural networks models; the first order equation; reviews. the Higuchi’s equation; the Peppas’ equation; In Table 1 the 153 works considered in this review are categorized the Hixon-Crowell’s equation; accordingly to the modeling approaches, whereas in Fig. 1.a a pie chart reporting the percentage of the modeling approach is shown. The ma- others. jority of works are based on simple fitting equations (∼46%), used to In Fig. 2 the percentage of use of these equations in the analyzed describe the drug release; the mechanistic models, which aim to works is reported. As it can be seen, the most used model (30%) is the mathematically describe the phenomenology of the HBSs, reach a bit one proposed by Peppas (comprising all the variants of the equation: less than one third of the total works (∼31%); statistical and neural see Section 2.2.2), followed by (19.4%) the Higuchi’s (which reduces to networks models are both at less than 10% of the total works and the be a particular case of the Peppas’ equation), the zero order model with reviews articles constitute ∼4% of the total database. Another inter- the 19% (constant release rate) and the first order drug release (15.6%). esting aspect is the number of publications per year, reported in The Hixon-Crowell’s equation constitute almost the 8% as well as all the Fig. 1.b. As it can be seen the total number of publications (black line other type of equations (mainly purely mathematical). with stars) regarding the “modeling of the drug release from hydrogels” The advantage of this type of models is surely the simplicity of is constantly increasing and the forecast, based on the last 5 years trend, utilization that has made them spread across the scientific community. confirms the growing trend. It is useful to have a deeper look at number The general form of these equations could be: of publications per year considering the works categorized on the Mt (t) modeling approach: the drug release fitting models (gray line with = kf(t) square), constitute the majority of the published work per year and M (2.1) have undergone a dramatic increase within the years 2000–2014, due where Mt is the total amount of drug released at time t , M is the total to the consolidations of the equations proposed in the late ’80 s and amount of drug released, k is a constant of the process representing the thanks to comprehensive reviews on the subject (i.e. Siepmann and HBSs’ characteristics (i.e. geometry, surface area etc…), f(t) is a generic Peppas (2001) Section 4). However since the 2014 the papers based on function of time. The release functionf(t) , depending on the model, can this approach have been diminishing (in favor of other approaches) and have “mechanistic” or empirical origin. In the following the most the forecasted trend predicts a further decrease. The number of pub- common release function will be derived, highlighting all the simpli- lications per year of mechanistic models shows that this approach (as fying assumptions that have to be made. considered in this review) has begun to be relevant in the early 2000s, where the peak could be attributed to the works of Siepmann et al. 2.1. Lumped parameters models (1999); Siepmann et al. (2000); Siepmann and Peppas (2000), and has maintained an almost constant trend until the 2014, where an uptrend The equations belonging to this category have been derived with was recorded and confirmed by the forecast. Statistical models are re- lumped element models, in which the delivery system is considered at ceiving much more interest in the last years, also thanks to the high uniform concentration. This can be valid for pure dissolving systems, in encouragement of the FDA to adopt the principle of quality by design which the internal concentration is constant, and in systems (device/ (QdB) in the new pharmaceutical product development (Yu et al., dissolution medium) where the mass transport limiting step is con- centrated in the outer dissolution medium. In these cases, the first being more physically sound than the latter, the drug release can be described 1 TITLE-ABS-KEY (((drug* AND release) OR (drug* AND deliver*)) AND with a lumped mass balance on the drug delivery device: (hydrogel* OR gel* OR crosslink* OR cross-link* OR tablet* OR hpmc) AND dW(t) (modeling OR modelling) AND NOT (pharmacokinetic* OR vivo OR rat OR pig dt = KS(csol d,sat csol d ) OR porcine OR bovine OR human OR patient OR clinic* OR cellular OR cells OR t=0 W(t) = W0 (2.2) cell OR bioavailability OR (molecular AND modelling) OR (molecular AND modeling) OR (fused AND deposition))) AND (EXCLUDE (DOCTYPE, "cp") OR where W(t) is the drug mass, K is the product of the mass transport EXCLUDE (DOCTYPE, "cr")) coefficient in the dissolving medium and the drug molar mass (Mw ), S is 177 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 Table 1 Articles dealing with “modeling drug release from hydrogels” in the time span 1980–2018, categorized accordingly to the modeling approach. Category Articles Drug release fitting models (Efentakis and Buckton, 1990; Ford et al., 1991; Skoug et al., 1993; Pather et al., 1998; Karasulu et al., 2000; Faisant et al., 2002; Buonocore et al., 2003; Khan and Craig, 2003; Gu et al., 2004; AlKhatib and Sakr, 2005; He et al., 2005; Moneghini et al., 2006; Reis et al., 2007; Blanco et al., 2008; Das and Senapati, 2008; Khan et al., 2008; Liu et al., 2008; Razzak et al., 2008; Zhang et al., 2008; Bera et al., 2009; Chakraborty et al., 2009; Ibn Razzak et al., 2009; Kabir et al., 2009; Kavitha et al., 2010; Odeku and Picker-Freyer, 2010; Shah and Shelat, 2010; Simi et al., 2010; Wadher et al., 2010; Hasanuzzaman et al., 2011; Tarirai et al., 2011; Ahad et al., 2012; Aslam et al., 2012; Ghosal et al., 2012; Jha et al., 2012; Kamlesh Jayantilal et al., 2012; Mohapatra et al., 2012; Naeem et al., 2012; Omprakash et al., 2012; Shahi et al., 2012; Deshpande et al., 2013; Ferrero et al., 2013; Hernawan et al., 2013; Mirzaei et al., 2013; Yin et al., 2013; Bhise et al., 2014; Dadkhah et al., 2014; Dini et al., 2014; Pan and Yang, 2014; Shah and Deshpande, 2014; Shah et al., 2014; Weiss et al., 2014; Ainurofiq and Choiri, 2015; Ganguly and Das, 2015; Gonullu et al., 2015; Letmanski et al., 2015; Maiti et al., 2015; Mauricio et al., 2015; Qazi et al., 2015; Tajrin and Mamun, 2015; Castillo-Miranda et al., 2016; Erum et al., 2016; Mabrouk et al., 2016; Mikac et al., 2016; Yadav et al., 2016; Aranaz et al., 2017; Cascone, 2017; Gaikwad et al., 2017; Heydari et al., 2017; Zandi et al., 2017; Gunda and Vijayalakshmi, 2018; Khalid et al., 2018; Khanum et al., 2018; Thapa and Jeong, 2018; Vidart et al., 2018) Mechanistic models (Peppas et al., 1980; Droin et al., 1985; Kou et al., 1990; Singh et al., 1994; Lu et al., 1998; Siepmann et al., 1999; Siepmann et al., 2000; Siepmann and Peppas, 2000; Faisant et al., 2002; Radu et al., 2002; Buonocore et al., 2003; Frenning and Strømme, 2003; Kiil and Dam-Johansen, 2003; Li et al., 2004; Siepmann et al., 2004; Yan et al., 2004; Frenning et al., 2005; Zhou et al., 2005; Agnihotri et al., 2006; Krishna Rao et al., 2006; Moneghini et al., 2006; Martinez et al., 2007; Reis et al., 2007; Abbasi et al., 2008a; Abbasi et al., 2008b; Saunders et al., 2008; Wu and Brazel, 2008; Siepmann et al., 2010; Wang et al., 2010; Zhu et al., 2010; Lamberti et al., 2011; Lamberti, 2012; Lamberti et al., 2012; Kan and Li, 2013; Kaunisto et al., 2013; Xu et al., 2013; Caccavo et al., 2015a, b; Lins et al., 2016; Salehi et al., 2016; Arroyo et al., 2017; Caccavo et al., 2017a; Caccavo et al., 2017c; Pareek et al., 2017; Siepmann et al., 2017; Caccavo et al., 2018; Ferreira et al., 2018) Statistical models (Tillotson and Sakr, 2004; Do et al., 2008; Kumari et al., 2009; Zu et al., 2011; Choi et al., 2012; Gaikwad et al., 2013; Wu et al., 2015; Iurian et al., 2016; Kumar et al., 2016; Shah et al., 2017; Sun et al., 2017; Barmpalexis et al., 2018; Farhadian et al., 2018; Yekpe et al., 2018; Zhang et al., 2018) Neural networks models (Ibric et al., 2002; Ibrić et al., 2003; Mendyk et al., 2006; Shao et al., 2007; Xie et al., 2008; Ivić et al., 2010; Chansanroj et al., 2011; Díaz-Rodríguez and Landin, 2012; Hezaveh et al., 2012; Petrović et al., 2012; Bagheri et al., 2014; Al-Zoubi et al., 2015; Yousefi and Razavi, 2017) Reviews (Siepmann and Peppas, 2001; Gabrielsson et al., 2002; Grassi and Grassi, 2005; Lin and Metters, 2006; Ganji and Vasheghani-Farahani, 2009; Zhang et al., 2012; Grassi and Grassi, 2014; Timmins et al., 2016) Fig. 1. a) The total number of scientific works (153) dealing with “modeling drug release from hydrogels” in the time span 1980–2018, categorized accordingly to the modeling approach. b) The number of publications per year (data smoothed with a central moving average on a time span of five years) dealing with “modeling drug release from hydrogels” in the time span 1980–2018, plus trends prediction (dashed lines), categorized accordingly to the modeling approach. the surface interested by the mass transport, csol d,sat is the drug saturation which is wrong for HBS as the others, but to the fact that a zero order concentration of the solution (it is supposed that the drug concentration release can be obtained in HBSs due to a combination of physical close to the device is constant and equal to the saturation condition) phenomena (i.e. diffusion and polymer relaxation) and therefore, in and csol d is the drug concentration in the dissolution medium. This mass light of that, its use could be considered legit. balance is also known as the Noyes–Whitney equation and it is ex- tensively used to describe dissolution phenomenon (with 2.1.1. The zero order equation K/Mw = Ddsol/h where Ddsol is the drug diffusion coefficient in solution The zero order equation, heavily used to describe drug release from and h the thickness of the boundary layer around the device). HBSs (Fig. 2) has the form: It has to be said that, a part from the zero order equation, the use of these models is conceptually wrong for hydrogels based systems since Mt (t) = kt they were derived to describe other phenomena (i.e. pure drug dis- M (2.3) solution) and, therefore, they are not related in any kind to the physics It can be easily derived from the lumped mass balances on the re- of the problem. The exception made for the zero order equation it is not lease device (Eq. (2.2)). In the special case of perfect sink condition, in related to its physical derivation done in the following paragraph, which the external concentration is negligible csol d = 0 (or more 178 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 whose solution gives Eq. (2.7). As it can be seen the Hixon-Crowell’s equation was derived for a completely different system (dissolution of pure drugs) and it constitutes only a particular solution (constant ex- ternal concentration and changing surface with a specific law). There- fore the Hixon-Crowell’s equation should not be used for HBSs. 2.2. Distributed parameters models The release kinetic from hydrogel based systems (Eq. (2.1)) is governed by the mass transport within the system itself, since the mass transport outside is order of magnitude higher and do not constitute a limiting step. In light of this, the release function f(t) can be derived from the solution of the continuity equation for the ith species (ne- glecting the convection transport) (Bird et al., 2007): ci = Ji t (2.9) where c i is the molar concentration of the ith species, is the nabla operator and J i is the molar diffusive flux. Eq. (2.9) holds true for any diffusive system, without any simplifying assumption. The choice of the form of the diffusive flux J i constitute a first as- Fig. 2. Percentage of use of the drug release fitting equations. sumption (1° assumption): the Fick’s law of binary diffusion that de- scribe the diffusive flux as J i = cD xi , where c is the total con- centration of the system, D is the diffusion coefficient and xi is the generally constant csol d = cost ), the integral of Eq. (2.2) gives: molar fraction of the ith species. This equation was derived to describe diffusion in dilute binary systems, where the total concentration c is W(t) W0 = KScsol d,sat t = k1t (2.4) constant (2° assumption) and the diffusion coefficient D is constant (3° Considering that the amount of drug in the device is related with the assumption), which, as it has been seen in the phenomelogy section, is amount of drug outside as: Mt (t) = W0 W(t) Eq. (2.3) can be easily very far from the real behavior. Eq. (2.9) can thus be particularized for obtained. the drug species (i = d) and rewritten as: cd 2c 2.1.2. The first order equation =D d (2.10) t The first order equation is frequently used to describe drug release from HBSs (Fig. 2), it has the form: This equation needs initial and boundary conditions to be solved. Usually the initial condition assumes uniform drug concentration (4° Mt (t) ln = kt assumption) and the boundary conditions consider constant drug sur- M (2.5) face concentration, usually zero in the perfect sink (5° assumption). It considers that the drug release is influenced by the external drug concentration, i.e. the perfect sink condition is not verified. Integrating 2.2.1. Higuchi’s equation eq. (2.2) considering that csol d = (W0 W(t))/V and supposing, for the Takeru Higuchi (1960) reviewed the state of art concerning the sake of simplicity, that csol d,sat = W0/V (the initial amount of drug in the percutaneous drug rate of penetration from ointments and creams. In volume V of dissolution medium is capable to reach the saturation case of the rate controlling step is inside the medicaments (not in the conditions): skin), both pure diffusion (“Absorption from Solutions”) and diffusion KS from drug suspension (“Absorption from Suspensions”) were con- W(t) = W0exp t V (2.6) sidered. While the first was later treated in detail by William I. Higuchi (Higuchi, 1962), the latter was better explain by T. Higuchi (Higuchi, Rearranging Eq. (2.6) considering that Mt (t) = W0 W(t) and 1961). Both approaches considered as limiting step the diffusion of drug M = W0 , the Eq. (2.5) can be easily be obtained. from the inside to the external medium, assuming a perfect sink con- dition on the external layer. 2.1.3. Hixon-Crowell’s equation In particular W.I. Higuchi considered the analytical solution of Eq. As it has been shown (Fig. 2), the Hixon-Crowell’s equation is (2.10) for a thin plane sheet (no edge effects), which is a one dimen- sometimes used to describe the drug release from HBS. It has the form: sional problem in rectangular coordinates (6° assumption), to derive W(t) 1 3 one of the most famous drug release fitting equation. Under these cir- =1 k1t cumstances Eq. (2.10) can be rewritten as: W0 (2.7) 2c It can be obtained from the solution of Eq. (2.2) in the special case cd =D d t x2 of perfect sink csold = 0 (or more generally constant cd = cost ), and sol t = 0, cd = cd0 considering that, according to Hixson and Crowell (1931), the surface x = ±l, cd = 0 (2.11) area changes proportionally to the drug dissolution with the function- 2 ality S = k sW(t) 3 , where k s is a constant containing shape and density The thin plane sheet has total thickness 2l and the problem is information. Eq. (2.2) can be rewritten as: symmetrical around the central plane of the sheet (the sheet can ex- dW(t) 2 change from both the surfaces). The analytical solution of this problem, = K1W(t) 3 dt in terms of fractional mass released is well known (Crank (1975), Eq. t=0 W(t) = W0 (2.8) 4.18): 179 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 layer in a planar geometry, he was able to relate the release, analo- gously to case of pure diffusion from solutions, to the square root of time. Also in this case the equation should be used with caution, since it is valid (considering all the other approximations) until the drug re- servoir (undissolved drug) is present. Therefore, while the “square-root” approximation could be used for already swollen systems (gel loaded) to describe the first part of the drug release, or for gel loaded with undissolved drug (that is readily dissolvable) until the drug reservoir is present, the Higuchi’s equation should not be used for swelling systems. 2.2.2. Peppas’ equations Nicholas A. Peppas has been one of the most prolific researcher in the field of hydrogels and HBSs for drug delivery. Among other inter- esting publications, his research group produced the most used equa- tion in the field: “the power law equation”, in different variations. In Korsmeyer and Peppas (1981); Langer and Peppas (1981); Korsmeyer et al. (1983) it was proposed a semi-empirical equation to describe the fractional drug release from thin plane sheets, much more flexible and adaptable to describe the experimental release data with respect to the Higuchi’s equation: Mt (t) = ktn M (2.15) Fig. 3. Fickian diffusion driven release kinetics from a thin plane sheet. D = 2 × 10 9m2 / s , l = 5mm. where n is an adjustable parameter ranging from 0.5 to 1 (for a plane 1 sheet geometry) and k = 4 ( ) , as demonstrated in the previous D (2l)2 2 paragraph. When n is 0.5, the Higuchi’s equation is retrieved, otherwise Mt (t ) 8 D (2n + 1)2 2 behavior different than the simple Fickian diffusion (which does not =1 exp t respect the assumptions on which the Higuchi’s equation is based) can M (2n + 1)2 2 (2l)2 (2.12) be described. When n is 1 the zero order release can be described, with n=0 In the short time limit the equation can be simplified as (Crank constant drug release rate (highly desirable to optimize the drug ad- (1975), eq. 4.20): sorption). 1 In Ritger and Peppas (1987a,b) the value of the coefficient n char- Mt (t ) Dt 2 1 n (2l) acteristic of the Fickian diffusion was obtained for geometries different =4 2 +2 ( 1)nierfc M 2l 2 n= 1 2 Dt (2.13) than the plane sheet. In particular it was analyzed the “infinite long” cylinder (with L ≫ D), which still verify the uniaxial mass transport (Eq. where ierfc (x ) is the integral error function Higuchi’s equation further (2.16)) and the sphere (Eq. (2.17)): simplify the solution dropping the summation part (7° assumption): Mt (t ) Dt 1 2 1 cd t = Dr 1 r (r ) cd r =4 = kt 2 t = 0, cd = cd0 M (2l)2 (2.14) r = ± a, cd = 0 (2.16) In Fig. 3 the solution of these equations for a model thin plane sheet is reported. As it can be seen both the simplified solution (Eqs. (2.13) and (2.14)) diverge sensibly form the exact solution (Eq. (2.12)) after cd t =D 1 r2 r (r ) 2 cd r t = 0, cd = cd0 the 60% of release. In Higuchi (1962) the author placed the safety limit for the use of Eq. (2.14) at max 30% of the total drug release. r = a , cd = 0 (2.17) Eq. (2.14) is of the most famous and cited equation to describe the where a is the radius of the cylinder and of the sphere. The exact and drug release from HBSs, it could be useful to recall the assumption on the short time solutions are known (Crank (1975), Eq. 5.23, 5.25 for Eq. which this equation is based: (2.16) and Crank (1975), Eq. 6.20, 6.22 for Eq. (2.17)): 1. Fickian diffusion; Mt (t) 4 2 =1 exp( D n t) 2. total concentration constant (dilute system); M n=1 a2 n2 (2.18) 3. constant diffusivity; 4. initial uniform drug concentration; Mt (t) =1 6 1 exp( D n2 2t/a2) 5. perfect sink conditions; M 2 n=1 n2 (2.19) 6. thin plane sheet geometry (with no edge effects); 1 3 7. short-time approximation (it can describe only the first 60% of the Mt (t) 4 Dt 2 Dt 1 Dt 2 = 1 + drug release); M 2 a2 a2 3 1 2 a2 (2.20) 8. constant geometry (non swelling systems) 1 Mt (t) Dt 2 1 na Dt =6 2 +2 ierfc 3 T. Higuchi (Higuchi, 1961), considered the drug diffusion from 2 M a n=1 Dt a2 (2.21) ointments with suspended drug (concentration greater than its solubi- lity), and assumed as limiting step the diffusion, not considering the where n are the roots of J0 (a n ) = 0 (J0 (x) is the Bessel function of the dissolution rate (instantaneous process). Considering a linear drug first kind of order zero). concentration profile from the internal drug reservoir to the external Ritger and Peppas (1987a) considered only the terms with time to 180 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 Fig. 4. a) Fickian diffusion from an infinite long cylinder (L > > D) and b) Fickian diffusion from a sphere. D = 2 × 10 9m2 / s , a = 5mm With the solid black line the exact solution in represented, with the dash blue line the short time approximation is reported, with the dot green line the term of the short time approximation function of the 0.5 power of time is reported and with the magenta dash-dot line the solution proposed by Ritger-Peppas. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the power of 0.5 to demonstrate that the Higuchi’s equation was not McNeill, 1984) could generate anomalous drug release. suitable to describe geometries different than thin plane sheet. To overcome this limitation it was proposed to keep the constant related to 3. Mechanistic models 1 the time to the power of 0.5 (k = 4(D / a2) 2 for cylinders and 1 k = 6(D / a ) 2 for spheres) obtained in the short time solutions and to 2 optimize the drug release modifying the exponent of the time. It was Analyzing the mechanistic works in can be noted that the HBSs are found as best fit the value of n = 0.45 for cylinders and n = 0.43 for mostly described with the “monophasic” approach (∼97%) and only spheres. The results of these approximations are show in Fig. 4. few works (Zhu et al., 2010; Xu et al., 2013) treat the problem under As it can be seen the solution proposed by Ritger-Peppas does not the conceptual assumption of “multiphasic” system. Within the mono- diverge as rapidly as the short time approximation simplified to one phasic category three different approaches can be individuated: term (in analogy to Higuchi’s equations) but neither describe very close the exact solution. It was also demonstrated that the mass transport the “single-component” approach (52%): takes into account only of from a tablets (independently on the aspect ratio) could be described the drug mass transport; with the sum of the short term approximation solution (superposition the “multi-components” approach with swelling approximation principle). (35%): consider the mass transport of, at least, water and drug, Ritger and Peppas (1987b) used the Eq. (2.15) to describe the drug describing the swelling with equations based on volumetric con- release from slightly swelling systems (equilibrium water content straints; ≤25%) with different geometries. The coefficient “n ” capable of de- the “multi-components” approach with mechanics (12.5%): con- scribing a constant drug release rate (zero order or Case II transport) sider the problem in its entirety striving to couple the mass transport was found: n = 1 for thin sheet layers, n = 0.89 for cylinders and with the mechanics (from which the swelling is retrieved). n = 0.85 for spheres. In Peppas and Sahlin (1989) the Fickian contribution, with a totally As it can be seen in Fig. 5, the single-component approach, heavily heuristic approach, was separated from the relaxation contribution: used at the beginning of the years two thousand have seen a sensible decrease, ranging now in a zone with less than one publication per year. Mt (t) The multi-components with the swelling approximation have seen a = k1tn + k2t2n M (2.22) peak in the year 2000, due to the “sequential layer model” proposed by Siepmann et al. (Siepmann et al., 1999; Siepmann et al., 2000; Being k1tn the Fickian contribution and k2t2n the relaxation con- Siepmann and Peppas, 2000), a slight decrease and a new uptrend tribution. It was postulated that, depending on the value of the fitting which is still present nowadays. The last approach, the most complete, parameters k1 and k2 (n are given from the previous works and depend is quite recent and it has become relevant in the last 3 years only on the geometry: 0.5 for plane sheets, 0.45 for long cylinders and (2016–2018). 0.43 for spheres) it was possible to estimate the relative contribution to the mass transport. Actually the use of Eq. (2.22), as well as the 2.15, to discriminate 3.1. The “single-component” approach between the types of mass transport, and therefore formulating phe- nomenological hypotheses, it is not advisable. Deviation from the In the single component approach the mass balance is performed on normal Fickian mass transport can depend on several variables and not the dissolved drug that diffuses within the system to reach the external only on the presence of diffusion coupled with the polymer relaxation. dissolution environment. In light of this, the release is governed by the The shape of the device, which is normally more complex than the continuity equation for the drug “d” species (neglecting the convection simple geometries analyzed, as well as the high degree of swelling (not transport) eventually considering a generation term (rd ) due to the drug accounted in these equations) and the initial drug profile (Graham and dissolution (Bird et al., 2007): 181 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 for water and drug, without considering the volume variations. Hariharan and Peppas (1993) modeled the mass transport of water, drug and ions in thin sheets (one-dimensional problem) considering as driving force the concentration gradients as well as the electrostatic potential and relating the swelling to the amount of water absorbed. Breakthrough works were published by Siepmann et al. (1999); Siepmann et al. (2000); Siepmann and Peppas (2000) (with almost 500 citations in the 2018) in which the “sequential layer” model was pro- posed. This model is based on the Fick’s law of diffusion with non- constant diffusion coefficients (of “Fujita type” (Fujita, 1961)) for the species water (i = w ) and drug (i = d ): ci t = (Di c i ) t = 0, ci = ci0 (3.3) With Dirichlet boundary conditions for both drug (perfect sink) and water (fixed water concentration). The model accounts also for the polymer erosion, with a release kinetic proportional to the exposed surface area (At ), considering a pseudo-steady state release of polymer: mp = mp0 k dissAt t (3.4) The model, solved numerically in 2D axisymmetric geometry, at each time step computes the concentration profiles of drug and water, by integration it calculates the amount (mass) of drug and water and, Fig. 5. The number of publications per year (data smoothed with a central from Eq. (3.4), the amount of polymer. Considering additive the specific moving average on a time span of five years) dealing with “mechanistic mod- volumes (ideally mixing), the total volume at each time step is calcu- eling of drug release from hydrogels” in the time span 1980–2018, plus trends lated and the new geometry is recalculated accordingly. This assumes prediction (dashed lines), categorized accordingly to the modeling approach. that the deformation is affine: i.e. a cylinder remain a cylinder within all the process. cd Buonocore et al. (2003) proposed a one dimensional model (thin = J d + rd t sheet) for water adsorption, drug release and swelling based on Fick’s t = 0, cd = cd0 (3.1) law for both water and drug. A peculiar boundary condition (BC) was used for the water concentration at the erosion front. It was considered With case specific boundary conditions. The majority (∼60%) of that the water concentration does not reach instantly the equilibrium works (Peppas et al., 1980; Lu et al., 1998; Faisant et al., 2002; Zhou concentration but it increases gradually depending on the relaxation of et al., 2005; Agnihotri et al., 2006; Krishna Rao et al., 2006; Martinez the macromolecules. It was used an empirical equation to calculate the et al., 2007; Siepmann et al., 2010; Wang et al., 2010; Kan and Li, 2013; water volume fraction (and therefore the water concentration used as Lins et al., 2016; Arroyo et al., 2017; Siepmann et al., 2017) considered BC) at the erosion front: very fast the dissolution compared to diffusion, therefore the drug is readily available to diffuse toward the external dissolution medium d w = (a1 w )[1 exp( (1 w ))] (rd 0 ), and constant diffusion coefficient, therefore obtaining from dt t = 0, =0 (3.5) the Eq. (3.1) the Eq. (2.10). As it has been shown previously this last w equation can be solved analytically, with solution equations depending where w is the water volume fraction and a1[=]1/s a constant to be on the geometry. determined. The results of this equation is a w (t) with sigmoidal shape, The remaining part (∼40%) used less limitative assumptions such useful to describe not immediate equilibrium condition at the interface. as a pure diffusive system with non-constant diffusion coefficient (Droin In Kiil and Dam-Johansen (2003) a one dimensional model (radial) et al., 1985; Kou et al., 1990) or diffusion (with non-constant diffu- of drug release from HBS was proposed. This model accounts for the sivity) influenced by dissolution (Singh et al., 1994; Frenning and presence of the erosion, diffusion and swelling fronts and their move- Strømme, 2003; Li et al., 2004; Siepmann et al., 2004; Yan et al., 2004; ments with empirical or mass constraint based equations. Frenning et al., 2005). The dissolution term, resembling the Noyes–- Grassi et al. (Grassi et al., 1998, 1999) first and Wu and Brazel Whitney equation, has the form: (2008) later, according to Camera-Roda and Sarti (1990), considered a molar flux different form the simple Fick’s law, to account for the in- rd = k( cd,sat c d) (3.2) fluence of the polymer relaxation on mass transport. where k is the first order dissolution constant, cd,sat is the equivalent ci drug saturation concentration in the hydrogel and cd is the field vari- t = (J di + Jri ) able. Eq. (2.10) with non-constant diffusion coefficient or Eq. (3.1) do t = 0, c i = c i0 (3.6) not have analytical solutions, and were solved numerically (mainly where J di is the Fickian flux and Jri is the polymer relaxation con- with finite difference schemes). tribution to the flux, with the form: Jri 3.2. The “multi-components” approach with swelling approximation Jri = Dir c i t (3.7) With this approach apart from the drug mass transport at least the where is the diffusion coefficient related to the polymer relaxation Dir water mass transport is modeled, considering the swelling/shrinkage of and is the polymer relaxation time. This equation is also known as the system function of the mass variations. Maxwell-Cattaneo equation, and it is used in the hyperbolic form of the Malley et al. (1987) solved the diffusion equations in thin sheets mass (and heat) transport equation to account for the finite speed of (one-dimensional problem) with Higuchi’s approximation (Eq. (2.14)) propagation of mass (and heat) signals within continuum fields. It was 182 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 proven to be able to describe the Fickian as well as the anomalous showing the year of publication, the total number of citations (which transport (Case II transport) (Camera-Roda and Sarti, 1990). Grassi can be seen as a measure of the spreading of the approach in the sci- et al. (1998, 1999) used the model (1-D) to show how the viscoelastic entific community), the type of geometry for which the modeling properties of the swelling matrix influence the features of the drug re- equations were particularized, the species considered, the form of the lease, resulting mainly in a delay in the drug release rather than in- mass flux, how the swelling was described, the possibility to describe ducing a shape modification of the release curve. Wu and Brazel (2008) erosion and the peculiarity of the work. The works of Siepmann et al. used a simplified version of Eq. (3.7), without the second term on the (1999); Siepmann et al. (2000); Siepmann and Peppas (2000) con- right-hand side, and were able to model (in a one-dimensional geo- stituted a breakthrough in the HBSs scientific community, presenting metry) anomalous water transport as well as the early time burst re- for the first time a model capable of describing diffusion and de- lease, calculating the time dependent dimensions on the basis of the formation in geometries different than the one-dimensional case. water presence in the rubbery region. However this model suffers from strong limiting assumptions such as In Lamberti et al. (2011), followed by (Lamberti, 2012; Lamberti affine deformations and not constrained polymer mass balance. The et al., 2012), it was used a diffusive flux (massic and not molar) for first produces a description completely different than the experimental water similar to the one of Wu and Brazel (2008), considering that the observations; the latter comes from the coupling of the erosion mac- characteristic time of diffusion is order of magnitude greater than the roscopic mass balance (for polymer) with the microscopic mass bal- polymer relaxation (in Eq. (3.7) 0 or analogously NDe. D 1). The ances (for water and drug) without physical based constraints (i.e. the deformation was still related to the water but not to its mass (volume sum of their mass/volume fraction equal to 1), leading eventually to integral of the concentration in the domain) but to its inlet flux and unphysical concentrations. All the following models worked to remove specifically to the “swelling flux” (Jri ): this allowed to describe non limiting assumptions, especially working on the form of the mass fluxes affine deformation in 3D geometries. It was postulated that the diffusive to deal with diffusion in relaxing polymer. The works of Buonocore flux J di is responsible for the inner core hydration, whereas the re- et al. (2003); Kiil and Dam-Johansen (2003); Wu and Brazel (2008) laxation (“swelling”) flux contributes to the hydration and volume pointed out some interesting aspects that could/should be considered in variation of the rubbery zone. With a local mass balance on the erosion modeling HBSs, such as a boundary condition different that a classic front (and not on the swelling front) the water mass fluxes were related Dirichlet for water, the possibility to track the erosion, diffusion and to local velocity of the domain boundaries. The erosion was considered swelling fronts or the use of more complex mass fluxes (such as the as an inward constant velocity (ker ). This, as shown in (Caccavo et al., Maxwell-Cattaneo equation) to account for anomalous transport. A part 2017a), by integration leads to a macroscopic polymer dissolution ki- from this last, the other considerations did not have a following in lit- netic of zero order proportional to the surface area: erature, due to newer experimental discoveries (i.e. confirmation of constant water concentration at the erosion front (Barba et al., 2009)) dmp dt = k 'er A (t) or due to difficulties related to their application in more complex t = 0, mp = mp0 (3.8) geometries. Lamberti et al. (2011) used the approach of Wu and Brazel (2008) to describe the water mass flux and removed the limiting as- where is a constant and A(t) is the surface area of the HBS. To ac- k'er sumption of affine deformations. Kaunisto et al. (2013), followed by count for the domain deformation along with the mass transport the Caccavo et al. (2015b), further refined the modeling approach reducing continuity equations for the species were coupled with an arbitrary the number of parameters required, constraining the polymer to the lagrangian eulerian (ALE) moving mesh. water and drug amount (mass fraction constraint) and improving the Kaunisto et al. (2013) analyzed the behavior of an HPMC matrix mass flux description. loaded with a poorly soluble drug, under the assumption of constant However it was also proved that with these models it is possible to density, coupling the polymer mass description with the transport describe the macroscopic water hydration and the drug release sa- equations for drug and water through the mass fraction constrain. Also tisfactorily well but completely missing the microscopic information in this case most of the work was done on the definition of the mass such as the concentration profiles and/or the deformation (Caccavo fluxes: the transport equations were based on a simplified version of the et al., 2017c). This pointed out two important aspects: experimentally it generalized Fick equation (Bird et al., 2007) where the driving force is is not sufficient to monitor the “macroscopic” drug and polymer release the species concentration gradient. Despite the elegant approach, all the and the water uptake to infer on the phenomenology of the system; the multicomponent interactions, except those with the solvent, were as- mass transport only based models (developed until now) tend to overfit, sumed to be zero and the multicomponent Fick diffusivities were in- possibly leading to erroneous interpretation of the phenomena. terpreted as “pseudo-binary”. For the water-polymer diffusivity a “Fu- jita-type” form was used. Like in Lamberti’s model, the swelling was described through an ALE moving mesh method but the swelling ve- 3.3. The “multi-components” approach with mechanics locity was derived from a polymer/solid drug mass balance on the erosion boundaries. As it has been discussed in the phenomenology section, the HBSs are Caccavo et al. (2015b), followed by (Caccavo et al., 2015a; Caccavo characterized by the influence of the mechanics on the mass transport et al., 2017a; Caccavo et al., 2017c), modeled the mass transport in HBS and vice versa, therefore a full mechanistic model cannot ignore this as a diffusive process using as mass flux the mixture average approx- peculiar aspect. In this area two types of category can be individuated: imation, derived from generalized Fick equation (Bird et al., 2007): the chemo-electro-mechanical models and the “fully-coupled” models. In Table 3 the models accounting for the mass transport and the me- xi M ji = i Di = Di i + Di i chanics in HBSs for drug delivery are reported. xi M (3.9) where xi is the molar fraction, i is the mass fraction, i the density of the species in the system ( i = i ) , M = ( i ( i /Mi ) ) 1 is the system 3.3.1. Chemo-electro-mechanical models average molar mass. The swelling of the system, as in Kaunisto’s and In this type of models, as it can be understood from the name, apart Lamberti’s models, was related to the mass fluxes, allowing to describe from the presence of water and drug, a certain relevance is given to ions non affine deformations by imposing a certain deformation velocity on (Saunders et al., 2008; Pareek et al., 2017). the erosion front. The model was proven to be able to describe de- The mass transport equation (continuity equation for the ith species) formations, polymer and drug release as well as hydration of HBSs. Eq. (2.9) is still valid, but the flux has to account for the electrical Summarizing, in Table 2 the most important models are reported, potential driving force for ionic species: 183  D. Caccavo Table 2 Most important “multi-components” models with swelling approximation. Year Work Citations* Geometry Species Flux Swelling Erosion Peculiarity 1998–1999 (Grassi et al., 1998, 1999) ∼50 1D 1. water Maxwell-Cattaneo molar flux with Related to the total No Maxwell-Cattaneo molar flux to account for 2. drug non-constant diffusivity amount of water polymer relaxation absorberd 1999 (Siepmann et al., 1999; ∼500 2D axisymmetric 1. water Fickian with non-constant Affine, related to the Macroscopic kinetic Diffusion and deformation (affine) in 3D Siepmann et al., 2000; 2. drug diffusivity total amount of geometry Siepmann and Peppas, 3. polymer (not species 2000) constrained) 2003 (Buonocore et al., 2003) ∼110 1D 1. water Fickian Related to the total No Buondary condition on non-constant water 2. drug amount of water concentration absorberd 2003 (Kiil and Dam-Johansen, ∼ 80 1D 1. water Fickian Related to the total No Tracking the evolution of the erosion, swelling 2003) 2. drug amount of water and diffusion fronts 184 absorberd 2008 (Wu and Brazel, 2008) ∼20 1D 1. water Simplified Maxwell–Cattaneo molar Related to the total No Simplified version of the Maxwell–Cattaneo 2. drug flux with non-constant diffusivity amount of water molar flux to account for polymer relaxation absorberd 2011 (Lamberti et al., 2011; ∼ 55 2D axisymmetric 1. water Simplified Maxwell–Cattaneo mass Related to the Microscopic (constant Non affine deformations Lamberti, 2012; Lamberti 2. drug flux with non-constant diffusivity “relaxation” water erosion velocity) et al., 2012), 3. polymer (not mass flux constrained) 2013 (Kaunisto et al., 2013) ∼20 2D axisymmetric 1. water Generalized Fick with non-constant Related to the polymer Microscopic (non- Generalized Fick for the mass fluxes and swelling 2. drug diffusivity mass flux constant erosion related to the polymer mass flux. Polymer mass 3. polymer velocity) fraction from the mass fraction constraint 2015 (Caccavo et al., 2015b, a; ∼ 60 2D axisymmetric 1. water Generalized Fick with mixture- Related to the polymer Microscopic (constant Generalized Fick with mixture-average Caccavo et al., 2017a; 2. drug average approximation and non- mass flux erosion velocity) approximation for the mass fluxes. Polymer mass Caccavo et al., 2017c) 3. polymer constant diffusivity fraction from the mass fraction constraint 4. excipients *(until October 2018). International Journal of Pharmaceutics 560 (2019) 175–190 D. Caccavo International Journal of Pharmaceutics 560 (2019) 175–190 Table 3 Models accounting for multicomponent mass transport coupled with mechanics. Year Work Citations* Geometry Species Flux Mechanics Erosion Type 2008 (Saunders et al., 2008) 8 2D 1. water Nernst–Planck Small No Chemo-electro- 2. drug deformations mechanical 3. polymer 4. ions 2017 (Pareek et al., 2017) 2 2D 1. water Nernst–Planck Small No Chemo-electro- 2. drug deformations mechanical 3. polymer 4. ions 2016 (Salehi et al., 2016) 4 1D 1. water Gradientof chemical potentials Small Yes Fully-coupled 2. drug deformations 3. polymer 2017 (Caccavo and Lamberti, 7 2D axisymmetric 1. water Gradientof chemical potentials Large No Fully-coupled 2017; Caccavo et al., 2018) 2. drug deformations 3. polymer 2018 (Ferreira et al., 2018) / 1D 1. water Gradient of concentrations (Fickian) plus Small No Fully-coupled 2. drug gradient of trace of stress tensor (non- deformations 3. polymer Fickian diffusion) *(until October 2018). Ji = Di c i Di F zi ci equation) is not influenced by the mechanics. Both the works (Saunders RT (3.10) et al., 2008; Pareek et al., 2017) used this model on an already (par- tially) swollen system, using constant diffusivities, analyzing the drug This is known as Nernst–Planck equation, it is used to describe the release and the swelling degree of the system, particularly focusing on motion of a charged chemical species in a fluid medium. It extends the ionic contribution (i.e. response to a pH variation). Fick's law of diffusion for the case where the diffusing particles are also moved with respect to the fluid by electrostatic forces. F is the Faraday constant, R is the gas constant, T is the temperature, zi is the valence 3.3.2. “Fully-coupled” models number of the ith species. It should be evident that for non-ionic species These types of models are nowadays the most advanced approach to (zi = 0 ) the Fick’s law is retrieved. Since the description of diffusion describe the HBSs’ behavior. In this review they have been called “fully- influenced by electrical fields is not the main aim of this review, the coupled” because they can describe the influence of the mass transport part on the calculation of the electrical potential field will not be on the mechanics and vice-versa, allowing to catch real behavior of treated here: readers interested in the subject could find useful the re- HBSs stimulated chemically and/or mechanically. ference (Li, 2010). As usual, the continuity equation for the ith species (eq. (2.9)) is still The structural field equation (force balance equation) of the system, valid whereas the molar flux has to present a more complex form to neglecting the inertial terms (quasi-static approximation) reads: account for the mechanical contribution. Salehi et al. (2016) showed that the in the “diffusional driving forces” (Bird et al. (2007) Eq. 24.1- =0 (3.11) 8), disregarding the presence of relevant body forces (such as electro- where is the stress field. The constitutive equation of are used to static potential field), the only driving force to the species diffusion is link the mechanics with the mass transport (and not vice versa, as it will the “concentration diffusion term”, whose driving force is the gradient be clear soon) and usually this models are solved in the assumption of of the chemical potential of the ith species, that cannot be simplified to small deformations, considering that is related to the deformation the gradient of the concentration (or mass/molar fraction) of the ith through an elastic modulus (E, Young’s modulus) and a Poisson’s ratio species (as in Fickian diffusion). The same results were obtained by ( ) (i.e. Hooke’s law for isotropic material) plus a volumetric term that Caccavo and Lamberti (2017); Caccavo et al. (2018) who derived the is demanded to describe the swelling/shrinkage ability: Posmotic. molar flux from the non-equilibrium thermodynamics and the dissipa- tion inequality (Gurtin et al., 2010). Therefore the diffusive flux (in ( ( u) Posmotic I) = 0 (3.12) absence of body forces) can be expressed as: where (u) is the stress tensor related to the deformation (u )